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X* Introduction 

The eventual aim of crop inventory studies is production 
estimation, not area or yield estimates alone. Production estimates can 
be made only at a level where area and yield strata intersect. The 
varianoe of the production estimates is dependent upon the means and 
variances of both area and yield in the stratum. Thus, it is important 
that the stratifications for area and yield estimation be coordinated, 
and that the levels for aggregation be selected so that acceptable 
variances are obtained. 

In the AgRISTARS Foreign Commodity Production Forecasting (FCPF) 
Project, estimates are to be made for corn and soybeans area and yield 
in the United States as well as in Brazil and Argentina. 

To make production estimates, NASA pr.. /ides area estimates based on 
analysis of remotely sensed data and the USDA provides yield estimated 
from a regression model > In order to obtain the most precise production 
estimates, the levels of estimation must be coordinated. Thus, a study 
to determine the precision at several possible levels of aggregation was 
proposed . 


2. Objectives 

The overall objective of this study was to determine the optimal 
level for combining area and yield estimates of corn and soybeans. 
Production estimates and their variances were computed for several 
levels of area and yield estimates, and the resulting estimates were 
compared. 


3,* Approach 

Iowa was selected to study the optimal level for combining area and 
yield estimates of corn and soybeans. This state was selected for study 
as it is included in the 1981 AgRISTARS pilot experiment. The year 
selected for evaluation (’’current year”) was 1978. 

The level at which aggregation of area and yield to obtain 
production should occur is dependent upon the technology being utilized 
for estimation. If, for example, area or yield estimates made at a 
given level are biased or unreliable, then aggregation at that level 
would most likely be undesirable regardless of any potential gains in 
precision. A change in the technology utilized for estimation, however, 
might produce reliable estimates at the same level and be a viable 
candidate for aggregation. This investigation assessed the optimal 
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level with respeot to the current technology. Current technology 
utilizes digital analysis of Landsat MSS data on sample segments to 
provide area estimates; regression models are developed from historical 
data and used with ourrent weather data to provide yield estimates. 
Several levels of obtaining both area and yield estimates were 
considered: county, refined stratum, reflned/split stratum, orop 
reporting district, and state. 

The model form and variables inoluded in the regression used by 
CCEA for yield estimation of corn and soybeans in Iowa were obtained. A 
weather data base with historical (at least 30 years) and "current year" 
weather data were needed for all the cooperative meteorological stations 
in Iowa. Historical and "current year" county area and yield estimates 
made by USDA/SRS in Iowa were acquired for the same time period. 

Coefficients for the regression equations were derived to predict 
yield using the historical weather and yield data at each of the levels 
of aggregation. A weather smoothing function was utilized to provide 
estimates of meteorological variables for the various strata studied. 
Using the 1978 weather data, "current year" yield estimates were made 
for corn and soybeans in Iowa. 

The yield estimate (Y) and its variance were computed based on the 
regression equations. The yield estimate was then aggregated to the 
state level using area weights. The aggregated yield variance was used 
to determine which stratification systems were candidates for precise 
estimation methods. 

For those levels of aggregation which appeared to be improvements 
over the currently used method, a further investigation into the effeots 
of using the current area estimation methodology was conducted. Within 
stratum variances for the area of the crops of interest were obtained. 
The production estimate (P) and its variance (V(P) ) were computed for 
all the candidate aggregations. Evaluations compared the variances with 
one another. 


3.1 Data Set Utilized 

For development of regression models for yield, a historical series 
of yield estimates and meteorological data were required. The USDA/SRS 
county level statistics for yield of corn and soybeans were obtained 
from the Iowa state of floe for 1932-78. The 1932-77 data were used in 
computing the regression coefficients, and the 1978 data were acquired 
for results comparison. 

Daily observations of temperature and precipitation for all the 
cooperative meteorological stations in the state of Iowa were purchased 
from the Iowa Geologioal Survey (1900-74) and some were supplied by 
another task (1975-78). 
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2*2 Levels of Aggregation 

During the Large Area Crop Inventory Experiment (LACIE), 
aggregation of area and yield estimates to production was done at the 
state level. Thus, this would be one level for investigation. 

P'or the state of Iowa, yield estimates will be made at the state 
level and one other level during the 1981 AgRISTARS pilot experiment. 
NASA/JSC requested that this level be the refined strata in the state 
(Figure i). The yield modeling group, however, thinking that these 
strata were too broad, suggested a subdivision of them (Figure 1). This 
subdivision will be referred to as the refined/spllt strata in this 
report. Both of these levels are being considered for evaluation. 

An additional level which seems to be natural to include is the 
crop reporting district level (Figure 2) a3 this has traditional], y been 
a standard unit for the reporting of agricultural statistics. Also, the 
county level is Included as the smallest possible unit using current 
yield estimation technology, as this is the smallest level for which 
historical yield estimates are available. 

Some characteristics of the strata are presented in Tables 1-3* 
Means and variability between counties within the strata are described. 


3.3 Meteorological Data Estimation 

In order to study the various levels of aggregation, yield 
estimates were needed at each of the levels. To make yield estimates 
using current technology, meteorological data were needed for each 
stratum. Not all counties contain weather stations; besides, weighting 
by nearby weather st^cons may provide a better estimate of the overall 
weather of a county than the use of one weather station alone (Figure 
3). 

For this reason, a weather smoothing routine was utilized. Wagner 
(1971) devised an objective analysis technique which incorporates a low 
pass filter and provides a good analysis in sparse data areas or with 
data containing significant noise. Furthermore, the characteristics of 
the applied filter function are easily calculated and the analysis 
technique is quite forgiving in terms of the sensitivity of choosing a 
filter function for a given data set. This technique was initially 
devised to remove high frequency fluctuations in the initial condition 
fields used for numerical weather forecasting. However, the consistency 
and speed of the technique make it a viable technique for our purposes. 

Odell (1975) compared ten techniques for interpolation for 
irregularly spaced sparse data: composite average, nearest neighbor, 
least squares linear regression, least squares convex hull, average 
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REFINED STRATA 



Figure 1, Maps of the refined strata developed at NASA/JSC (top) and 
the refined/split strata as subdivided for the yield modeling effort 
( bottom) . 
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TabJ 1. Some characteristics of the refined strata. Means and 
variability are described for oorn and soybeans proportions and 
yields. 


DESCRIPTION OF PROPORTIONS FOR REFINED STRATA 


Stratum 


CORN 



SOYBEANS 


No. of 
Counties 

Mean 

Standard 

Deviation 

c.v. 

Mean 

Standard 

Deviation 

c.v. 

14 

37.0 

4.3 

11.6 

16.7 

3.8 

22.6 

13 

24 

37.7 

9.0 

35.4 

25.5 

9.0 

35.4 

45 

25 

25.7 

8.7 

33.7 

13.3 

5.0 

37.3 

41 


DESCRIPTION OF HELDS FOR REFINED STRATA 



CORN 



SOYBEANS 





Standard 



Standard 


No. of 

Stratum 

Mean 

Deviation 

C.V. 

Mean 

Deviation 

C.V. 

Counties 

14 

85.1 

23.4 

27.5 

32.3 

3.9 

12.1 

13 

24 

99.6 

15.7 

15.8 

32.6 

4.0 

12.3 

45 

25 

94.7 

18,2 

19.2 

31.0 

4.6 

14.8 

41 
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Table 2. Some characteristics of the refined/spllt strata. Means and 
variability are described for oorn and soybeans proportions and 
yields. 


DESCRIPTION OF PROPORTIONS FOR REFINED/SPLIT STRATA 


Stratum 


CORN 



SOYBEANS 


No. of 
Counties 

Mean 

Standard 

Deviation 

c.v. 

Mean 

Standard 

Deviation 

C.V. 

14A 

39.3 

2.9 

7.5 

14.1 

2,6 

18.7 

6 

14B 

35.0 

M.M 

12.6 

19.0 

3.1 

16,4 

7 

24A 

39.7 

2.0 

5.1 

29.7 

6,6 

22.1 

13 

2MB 

39.7 

3.2 

8.1 

29.3 

5.4 

18.4 

20 

2MC 

32.3 

5.2 

16.1 

14.6 

7.1 

48.9 

12 

25k 

21.3 

8,8 

Ml ,1 

12.4 

3.5 

28.4 

18 

25B 

22.2 

5.7 

25.8 

0,9 

0.2 

20.2 

3 

25C 

30.2 

C.7 

22.2 

16,0 

2.9 

18.0 
if i 0* Hi *m tm • 

20 


DESCRIPTION OF YIELDS FOR REFINED/SPLIT STRATA 


Stratum 


CORN 



SOYBEANS 


No. of 
Counties 

Mean 

Standard 

Deviation 

C.V. 

Mean 

Standard 

Deviation 

C.V. 

14A 

81.7 

22.9 

28.0 

32.2 

4.8 

14.9 

6 

1 MB 

88.0 

23.7 

26.9 

32.4 

3.1 

9.6 

7 

24A 

97.5 

18.5 

19.0 

33.5 

4.1 

12.2 

13 

24B 

103.6 

14.1 

13.6 

33.0 

3.8 

11.5 

20 

24C 

95.4 

13.2 

13.8 

30.8 

3.9 

12.7 

12 

25A 

85.6 

20.4 

23.8 

29.1 

4.3 

14.8 

18 

25B 

102.4 

12.5 

12.2 

32.9 

4.2 

12.8 

20 

25C 

98.6 

8.4 

8.5 

30.0 

2.8 

9.3 

3 
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Table 3. tfome characteristics of the orop reporting districts, Me. na 
and variability are deaoribed for oorn ana soybeans proportions and 
yields. 


DESCRIPTION OF PROPORTIONS FOR CROP REPORTING DISTRICTS 




CORN 


• iO ■**■»*•*• 

SOYBEANS 



Stratum 

Mean 

Standard 

Deviation 

C.V. 

Mean 

Standard 

Deviation 

C.V. 

No. Of 
Counties 

North West 

39.5 

2.3 

5.7 

26.8 

8.5 

31.7 

12 

North Central 

39.9 

2.1 

5.2 

30.7 

5.4 

17.5 

11 

North East 

30,2 

7.1 

23.5 

11.7 

8.9 

75.7 

11 

West Central 

38.3 

4.7 

12.4 

19.3 

8.3 

43.1 

12 

Central 

37.6 

4,9 

13.1 

25.1 

6.5 

25.9 

12 

Eaat Central 

33.1 

5.3 

16.0 

13.5 

3.4 

40.1 

10 

South West 

29.7 

f J,) 

19.6 

17.5 

4.0 

22.8 

9 

South Central 

16.3 

4.S 

29.3 

10,6 

2.3 

22.0 

11 

South East 

27.0 

6.9 

25.7 

16.5 

3.1 

18.7 

11 


DESCRIPTION OF YIELDS FOR CROP REPORTING DISTRICTS 


Stratum 


CORN 



SOYBEANS 


No. of 
Counties 

Mean 

Standard 

Deviation 

C.V. 

Mean 

Standard 

Deviation 

C.V. 

North West 

93.1 

21.0 

22.6 

33.8 

4.3 

12.7 

12 

North Central 

99.8 

14.1 

14.1 

32.2 

3.5 

10.9 

11 

North East 

95.9 

12.7 

13.2 

29.8 

3.5 

11.7 

11 

West Central 

89.9 

21.0 

23.4 

31.8 

3.9 

12.3 

12 

Central 

106.8 

12.8 

12.0 

33.9 

3.8 

11.2 

12 

East Central 

100.5 

12.7 

12.6 

33.8 

3.6 

10.7 

10 

South West 

85.7 

25.1 

29.3 

31.2 

3.4 ' 

10.9 

9 

South Central 

86.0 

18. 7 

21.7 

28.2 

4.7 

16.7 

11 

South East 

101.9 

12.8 

12.6 

31.7 

4.4 

13.9 

11 


X 


X 



igure 3* An example of a situation when weighting by weather stations 
in adjacent counties may be beneficial in providing good estimator, of 
weather for county k. Each X represents a meteorological station. 
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linkage, average linkage with directional correlation, Wagner* a 
objective analyaia, modified linkage, and modified least squares. These 
techniques were tested in terms of their ability to Interpolate five 
years of wheat yield data aoross the state of North Dakota (45 data 
points) based on seven stations of wheat yield data* The weighted 
linear regression technique appeared to bo the best technique with the 
objective analysis, least squares linear regression, and the modified 
average linkage coming In close behind. However, the weighted linear 
regression is computationally time consuming, the least squares linear 
regression is not well behaved on the boundaries, and the modified 
linkage does not m fleet directional trends in the data. The objective 
analysis approach provides a smooth well-behaved surface and is 
computationally fas . Its major deficiency is that the original data 
points are not fit exactly. However, if noise exists in the input data, 
this oan be advantageous, And, use of the cooperative meteorological 
station data makes this a reasonable assumption. 

Integration of data fields (raster form) produoed by the objective 
analysis routine is sometimes required in order to obtain averages of 
meteorological (or other) data over some polygonal area. In order to 
accomplish this, the subroutines of Rios (1979) were utilized. A driver 
program was written to enable averages, mean square errors, and 
variances to be calculated for polygonal areas with 39 or fewer 
vertices. The polygon may contain both convex and concave features. 
This capability enables averages for a farmer's field, or an entire 
political subdivision or stratum to be calculated. 

The general procedure utilized by the objective analysis technique 
is Illustrated by Figure 4. A grid of a user-selected density is placed 
over the area of interest. Then the available met station data are used 
to specify the values at the nearest grid intersection points. The 
objective analysis procedure then uses gradient and Laplaolan weights to 
speoify the values at all grid intersections (Wagner, 1971). Finally, 
an estimate of the smoothed variable can be made over any polygon of 
Interest by averaging over the grid points within that polygon. 

The objective analysis technique was found to perform well in 
interpolating maximum temperature, minimum temperature, and 
precipitation on both a monthly and a daily basis for a case study in 
May 1977 in Oklahoma (Pitts, 1980). 

Based upon the favorable results obtained by other investigators, 
the FORTRAN coded programs for objective analysis were obtained from Dr. 
David E. Pitts of NASA/JSC. The programs were modified to fit our 
specific needs; the resulting listings are presented in Appendix A. 

A meteorological data smoothing experiment was conducted to 
determine how the objective function should be utilized. One month of 
daily data (June 1974) for all met stations in Iowa was used in the 
study. There were several factors in the experiment: grid size (25 x 



n 





Nearest Grid intersection 



Objective Analysis 
Specifies All Grid Points 



Within Polygons of Interest 


Figure 4. Schematic diagram of the steps in the meteorological data 
smoothing routine used to obtain meteorological estimates for 
polygons of interest in Iowa, 
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25, 32 x 32, 64 x 64), level of smoothing (dally vs. monthly), gradient 
weight (1,10), and Laplaoian weight (1,10), The results were evaluated 
by examining the mean square error of fit to station data and the 
maximum change in specified valuoe. 

The first observation from this experiment was that using gradient 
and Laplaoian weights of 10 caused too much change in the specified 
values. A difference of up to about one lnoh of precipitation was 
observed. Thus, the remainder of the experiment was analyzed using 
weights of one only. 

The maximum absolute deviation from specified values was examined 
for the three grid sizes (Table 4), The 64 x 64 grid provided estimates 
much oloser to the specified values than the other two grid sizes. The 
root mean square error was examined for daily vs. monthly averaging 
(Table 5). It was found that averaging met data to monthly values and 
then smoothing the monthly averages performed significantly better than 
smoothing dally values and then averaging the smoothed values to obtain 
a monthly estimate. 

The parameters selected for use in our study were: grid size 64 x 
64 over Iowa, gradient and Laplaoian weights of 1.0, and smoothing of 
monthly average values. 


3.4 Yield and Yield Variance Estimation 

Estimates of yield at all the levels of aggregation were required 
for this study. To do this, the variables used in the CCEA state level 
model were utilized (Table 6). Regression coefficients were developed 
for each set of strata utilizing 1931-77 meteorological data and 1932-77 
USDA/ESCS estimates of county level yields. Data from 1970 were not 
used in the derivation of regression coefficients for corn due to the 
corn blight which occurred during that year. The meteorological data 
inputs were daily reports of minimum temperature, maximum temperature, 
and precipitation for the strata as computed by the Wagner variational 
analysis technique. 

The programs used to compute the yield and the yield variances were 
written in SAS. A sample program for each crop is given in Appendix B, 

The variance of the yield estimates was computed from the 
regressions by: 

V(?) - a 2 . (1 + x’ (xMx'x)' 1 *). 
y x 

Using this formula, a variance was computed for each of the strata in 
each of the candidate stratification systems. 


Table Some results from the meteorological data smoothing 

experiment. The table shows daily maximum absolute deviations of 
smoothed values from the specified station values. 


Weather 

Variable 


Grid Size 


25 x 25 

32 x 32 

64 x 64 

Maximum Temperature 

2.93 

2.45 

0.77 

Minimum Temperature 

2.08 

1.39 

0.63 

Precipitation 

0.06 

0.04 

0.01 


Table 5. Some results from the meteorological data smoothing 
experiment. The table shows the root mean square error of smoothed 
values from the specified station values. 


Weather 

Variable 

Grid 

Size 

RMS 

Error 

Daily Smooth 

Monthly Smooth 

Temperature 

32 x 32 

4.88 

0.52 


64 x 64 

NA* 

0.17 

Precipitation 

32 x 32 

6.67 

0.92 


64 x 64 

NA 

0.17 


* NA - Not Available 
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Table 6. Model variables for the regressions predicting yield of corn 
and soybeans in Iowa. 


Corn 

Soybeans 

Linear trend 1941 >0 

Linear trend 1932-74 

Linear trend 1961-72 

Cumulative preoipitation 
October - April DFN 

May temperature x 

precipitation interaction 

May temperature x 

preoipitation interaction 

June temperature x 

precipitation interaction 

June temperature DFN* 

2 

June temperature (DFN) 

July precipitation DFN 

July precipitation DFN 

July temperature DFT 

July temperature DFT 

August precipitation DFN 

2 

July temperature (DFT) 

p 

August precipitation (DFN) 

August temperature DFT 

August temperature DFT 


* DFN = Departure From Normal 
DFT = Departure From Trend 



15 


To compare the preolalon of the several levels, a yield estimate 
aggregated to the state level was computed. The variance of the 
aggregated estimates must account for both the variance and covariances 
of the estimates. Thus, if 

$L " l w i 

then its variance can be computed by: 

V(Y l ) « ! W l VCYi) + 2 i E j Cov (Yi^j) WiWj 

where Yl is the aggregated state yield estimate, w* is the area weight 
for stratum i, and is the regression yield estimate for stratum i. 


3.5 Area and Area Variance Estimation 

*■ 

The area estimates used in the study were the 1978 final area 
harvested estimates made by the USDA/SRS for the Iowa counties. The 
variance of these estimates is not computed due to the complex 
estimation methodology employed. More Importantly, the variance of the 
USDA estimates would most likely differ from that obtained utilizing the 
AgRISTARS "current technology” of estimating crop areas based on 
analysis of Landsat data over sample segments. 

The variance of the area estimates was computed using the methods 
described by Chhikara and Perry (1980). The number and distribution of 
agricultural segments in a region were obtained from NASA. For regions 
without Landsat imagery, the unobserved potential segments were assigned 
the same distribution of percent agricultural as the observed segments 
in that county. 

Two methods (Chhikara and Perry, 1980) were available to fit the 
model : 

o 2 * Ax B , 
x 

one method based on field sizes and the other a pixel-based method. 
Both methods were used for comparison and verification. 

In the field size model: 

a x a * I ? 

was computed. The field size estimates were obtained from Pitts (1980) 
data base. Counties were assigned an average field size equal to that 
of any sample segments within the county. Counties without segments 

were assigned the field size of a county with similar farm size in 
geographic proximity. 
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The second computational method, based on pixels, used the 
equation: 

al o ■ +c»2^ +a 3 ( 0 * 3^82 - i^ + P ). 

The programs which were utilized to estimate area varianoes are 
presented in Appendix C. 


3,6 Aggregation Methodology 

For each of the strata in each stratification system, the 
production was oomputed as: 


AAA 

P i * A i x Y i 

where P* , A* , and Y* are the estimates of production, area, and yield, 
respectively, in stratum i. 

The state-level aggregated production estimate was computed using 
the, methods of stratified estimation presented in Cochran (1963)* 

For each stratum, the variance of production was computed as: 

V(p i ) - V(& ± x Y ± ) - VCA^ + + VCAiWY-t) 

For the aggregated state estimate, 

V(P L ) - V(S P^ - EVO^) + 2^ Cov (P it Pj). 

4. Results and Discussion 


4.1 Regression Analysis 

A selection of results from the regression analyses are illustrated 
in Figures 5-9. At the county level, the correlations between values 
predicted by the regression and USDA observed values had a substantial 
range. Linn County, with r-square values of 0.93 and 0.92 for corn and 
soybeans, respectively, is fairly representative of a high correlation 
situation. Lyon County (r-squares of 0.78 and 0.76) is representative 
of a lower correlation. 

In estimating the yield of corn and soybeans at the crop reporting 
district level, the correlations between estimated and observed values 
did not have such a large range as for counties. This is due (at least 
in part) to the smoothing effects achieved by the consideration of a 
larger geographic region. The results for the state are illustrated 
also. 


I 
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IOWA CORN MODEL 

UNN COUNTY 
R-SOUARE « .030 



LEGEND CODE 4 * * USDA ESTIMATE « ° <> PREDICTED 


IOWA SOYBEAN MODEL 

UNN COUNTY 
R-SOUARE > .018 



YEAR 

LEGEND CODE * » * usd* ESTIMATE o o o PREDICTED 


Figure 5. Comparison of corn and soybean yields predicted by the 
regression equations with USDA/SRS estimates for Linn County, Iowa. 
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IOWA COHN MODEL 

LYON COUNTY 
R-SOUARE - .777 



30 33 36 30 42 45 48 51 54 57 60 63 66 60 72 75 70 


YEAR 

LEGEND CODE * * * USDA ESTIMATE » • * PREDICTED 

IOWA SOYBEAN MODEL 


LYON COUNTY 
R-SOUARE - .761 



30 33 36 30 42 45 40 51 54 57 60 63 66 60 72 75 70 

YEAR 

LEGEND. CODE ♦ * • USDA ESTIflATE • ♦ * PREDICTED 


Figure 6. Comparison of corn and soybean yields predicted by the 
regression equations with USDA/SRS estimates for Lyon County, Iowa. 
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IOWA CORN MODEL 

NORTH WiST CROP REPORTINC DISTRICT 
R-SQUARE > 866 



30 S3 36 39 42 45 48 51 54 57 60 63 66 66 72 75 78 


YEAR 

LECENDt CODE • * * USDA ESTIHATE » « * PREDICTED 
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Figure 7. Comparison of corn and soybean yields predicted by the f 

regression equations with USDA/SRS estimates for the North West Crop fi 

Reporting District in Iowa. i 
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Figure 8. Comparison of corn and soybean yields predicted by the 
regression equations with USDA/SRS estimates for the East Central 
Crop Reporting District in Iowa. 
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4.2 Variance of Yield Estimates 

The variance of the yield estimates was computed for each stratum 
in each of the stratification systems from the regression equations. 
The aggregate of these results to the state level is shown in Table 7 
and illustrated in Figure 10. 

For both crops, the state level had the largest variance, and the 
county level had the smallest variance. The refined strata had somewhat 
larger variances than those associated with the reflned/spllt strata. 
The variances for the CRD and the refined/split strata are about the 
same. 


4.3 Variance of Production Estimates 

The variance of the production estimates, computed based on the 
proceeding results, is shown in Table 8 and is illustrated in Figures 11 
and 12. There were only small differences between the values computed 
using the field size and pixel size methods for computing area 
variances. Thus, the same discussion applies no matter which method is 
selected. 

For corn, the difference in the standard deviations among levels is 
not great. The differences for soybeans are quite apparent, however, 
with the standard deviation at the state level being more than 50 t 
greater than for the other levels. For both crops, the refined strata 
had the smallest standard deviations. 

This result is somewhat surprising since the aggregated yields had 
shown this method to be of slightly lower precision. What is probably 
being illustrated, however, is the precision gained by having fewer 
strata. Ther® are only three refined strata compared with eight 
refined/split strata and nine crop reporting districts. Due to the 
strata correlations, relatively more precision was obtained with the 
fewer refined strata. 


5. Summary and Conclusions 

Aggregation of area and yield to production at the state level was 
the least precise of the methods examined. The crop reporting district, 
refined strata, and refined/split strata had similar levels of precision 
of production estimates. 

In examining the variance of yield estimates, the aggregated 
results from estimation at a county level showed a high precision. 
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ESTIMATED 1978 CORN YIELD FOR IOWA 



ESTIMATED 1978 SOYBEAN YIELD FOR IOWA 



Figure 10. The estimated yield of corn and soybeans at the state level 
for each of the stratification systems. The shaded area is the 
estimated yield plus and minus one standard deviation. 
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Table 8. Standard deviation of state-level production estimates 
made at several levels of aggregation. Units are In thousands 
of bushels. 


Level of 
Aggregation 

Field Size Method 

Pixel Size Method 

Com 

Soybeans 

Corn 

Soybeans 

CRD 

72.2 

14.8 

72.2 

14.8 

Refined 

63,8 

13.8 

63. 3 

13.6 

Refined/Split 

71.1 

14.2 

71 „0 

14.1 

State 

71.5 

22.5 

69.2 

22.2 
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ESTIMATED 1978 CORN PRODUCTION FOR IOWA 

FIELD SIZE VARIANCE ESTIMATION METHOD 



ESTIMATED 1978 SOYBEAN PRODUCTION FOR IOWA 

FIELD SIZE VARIANCE ESTIMATION METHOD 



Figure 11. Estimated corn and soybean production for Iowa using the 
field size variance estimation method. Shaded area is estimated 
production plus and minus one standard deviation. 
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ESTIMATED 1978 CORN PRODUCTION FOR IOWA 

PIXEL SIZE VARIANCE ESTIMATION METHOD 



ESTIMATED 1978 SOYBEAN PRODUCTION FOR IOWA 

PIXEL SIZE VARIANCE ESTIMATION METHOD 


300000 


200000 


1 00000 



Figure 12. Estimated corn and soybean production for Iowa using the 
pixel size variance estimation method. Shaded area is estimated 
production plus and minus one standard deviation. 
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Unfortunately, we did not have a mechanism for estimation of area 
variances at the county level. This level for estimation should be 
considered, however, because of its high precision for yield. 

The results of this study provide a first step in determination of 
the optimum level for combining area and yield estimates to obtain 
production. Two aspects not considered in this study should be part of 
a further analysis: (1) bias introduced using area or yield estimation 
at any of the levels, (2) the cost of computation of area and yield 
estimates at the varying levels, and (3) differences in the optimal 
level due to the crop of interest. 

The bias of estimates was considered in this study only to the 
extent that it did not appear that any of the estimates were biased with 
respect to the estimates made using any other level of estimation. It 
may be, however, that either area estimates, yield estimates, or both 
may have a bias when estimated at one of the levels. The area estimates 
are currently made at a refined stratum level; no information is 
available on the potential bias introduced by estimating areas on any 
smaller geographic region. The yield estimates are now made generally 
at the state level. Biases may be introduced due to the density of 
weather stations available for estimating the parameters of the 
regression equation. A technique such as was utilized in this study may 
be one possible solution to this problem. However, it is possible that 
the resulting yields should be smoothed rather than the input 
meteorological data since the relationship between the input data and 
predicted values is not linear in the input variables. 

The costs of computing the area and yield components must also be 
considered before a final recommendation can be made. The basis for the 
decision will consist of consideration of the variances and standard 
deviations computed as a part of this investigation coupled with cost 
information for computation of area and yield estimates at each of the 
potential levels of aggregation. The analysis can be carried out based 
on sample survey design theory such as described by Cochran (1963). 

The results for corn and soybeans were substantially different, 
with the level at which corn is aggregated making less difference than 
the level at which soybeans are aggregated. Thus, additional crops of 
interest such as small grains should be examined. 

In summary, the results of this study indicate that aggregations 
should be performed at a level below the entire state. Selection of the 
most appropriate level, however, requires further study of bias, cost, 
and crop-dependent differences. 
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FORTRAN programs used to carry out the meteorological 
routine. 



ooo ooonoooooooooooooooooooooooooooo 


31 


SMOOTH FORTRAN 

WRITTEN BY DAVE PITTS (J.S.C.) 

MODIFIED BY KEVIN MCCULLEN (L.A.R.S.) 06/09/80 

SMOOTH FORTRAN IS THE INITIAL CALLING PROGRAM FOR A WEATHER DATA 
SMOOTHING FUNCTION 


VARIABLES 


MAXLAT 

MINLAT 

MAXLNG 

MINLNG 

ISIZE 

JSIZE 

L 

IDEBUG 
IBUF(2) 
I POLY 


MAXIMUM LATITUDE ON MAP 

MINIMUM LATITUDE ON MAP 

MAXIMUM LONGITUDE .ON MAP 

MINIMUM LONGITUDE ON MAP 

NUMBER OF GRID IN NORTH-SOUTH DIRECTION 

NUMBER OF GRID IN EAST-WEST DIRECTION 

NUMBER OF PAIRS OF VERTICES OF POLYGON 

EQUALS L IF EXTRA PRINTOUT IS NEEDED FOR DEBUGGING 

NUMBER OF PSEUDOZONE OR FIELD DESCRIBED BY POLYGON 

NUMBER OF PbOYGONS TO BE PLACED OVER THE FUNCTION U 


ISIZE, MAXLAT, AND MINLAT MUST BE ADJUSTED SO THAT ISCALE IS AN 
INTEGER. 

JSIZE, MAXLNG, AND MINLNG MUST BE ADJUSTED SO THAT JSCALE IS AN 
INTEGER. 


IMPLICIT INTEGER**! (I-N), REAL* 8 (A-H, O-Z) 

REAL*8 XLAT(500) ,XL0NG(500) ,TMAX(500) ,TMIN(500) ,PREC(500) 

REAL* 8 P(64,64),U(64,64),DIFF(64,64) 

REAL*8 MAXLAT, MINLAT, MAXLNG, MINLNG 
INTEGER *4 IBUF(80) ,1X6(512) 

INTEGER *4 NTIMES 

5 READ (5,101 ,END=9200) MAXLAT, MINLAT, MAXLNG, MINLNG, ISIZE, JSIZE, K, 
+AA , ALF2 , ALF4 , ERR , MAXPAS 

101 F0RMAT(F9.3,$F10.3,3I5,/,4F10.3,I5) 

ISIZE = 64 
JSIZE = 64 
NTIMES = 64 

ISCALEs ( FLOAT ( ISIZE- 1 ) ) / (MAXLAT-MINLAT ) 

JSCALEs ( FLOAT ( JSIZE- 1 ) ) / (MAXLNG-MINLNG ) 


SM000010 
zSM000020 
SM000030 
SM000040 
SM000050 
SM000060 
SM000070 
SM000080 
SM000090 
SM000100 
sSMOOOl 10 
SM000120 
SM000130 
SM000140 
SM000150 
SM000160 
SM000170 
SM000180 
SM000190 
SM000200 
SM000210 
SM000220 
SM000230 
SM000240 
SM000250 
SM000260 
SM000270 
SM000280 
SM000290 
SM000300 
=SM000310 
SM000320 
SM000330 
SM000340 
SM000350 
SM000360 
SM000370 
SM000380 
SM000390 
:=SM000400 
SM000410 
SM000420 
SM000430 
SM000440 
SM000450 
SM000460 
SM000470 
SM000480 
SM000490 
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SM000500 

:n=s READ IN LOW DENSITY MAP DATA ================================r====SM000510 

SM000520 

CALL IN ( XLAT , XLONG ,TMAX, TMIN , PREC , L , MAXLAT , MINLAT , MAXLNG t MINLNG ) SM000530 

DO 104 I=1,NTIMES SM000540 

DO 10H Jsl ,NTIMES SM000550 

104 P(I,J)sO.O SM000560 

DO 102 M=1,L SM000570 

XI=(MAXLAT-XLAT(M) )»ISCALE+1 SM000580 

XJ = ( MAXLNG -XLONG (M ) )«JSCALE+1 SM000590 

IsXI SM000600 

IF (XI-FL0AT(I).GE.0.5) IsI+1 SM000610 

J=XJ SM000620 

IF (XJ-FLOAT(J) .GE.0.5) JsJ+1 SM000630 

IF (K.EQ.2) P(I, J)sTMAX(M) SM000640 

IF (K.EQ.O) P(I,J)sTMIN(M) SM000650 

102 IF (K.EQ.1) P(I, J)=PREC(M) SM000660 

SM000670 

;==s WRITE FIELD, DO OBJECTIVE ANALYSIS, PRINTOUT CONTOURED RESULTS ===SM000680 

SM000690 

WRITE (6, 3) SM000700 

3 FORMAT(1H1,'64 X 64 GRID*) SM000710 

WRITE (6,2) (<P(I,J),Js1,l6),Is1,NTIMES) SM000720 

WRITE (6,2) ((P(I,J),J*17,NTIMES),I=1,NTIMES) SM000730 

2 FORMAT (32(16(1X*F7.3),/)) SM000740 

CALL ANALdSIZE , JSIZE , AA ,ALF2 , ALF4 ,MAXPAS ,ERR ,P ,U ) SM000750 

WRITE (6, 3) SM000760 

WRITE (6,2) ((U(I,J),J*1,l6),Is1,NTIMES) SM000770 

WRITE (6,2) ((U(I,J),J=17,NTIMES),I=1,NTIMES) SM000780 

DO 9999 1:1,64 SM000790 

DO 9999 Jsl ,64 3M000800 

DIFF(I, J)sP(I, J)-U(I, J) SM000810 

999 IF (P(I,J).EQ.0.0) DIFF(I,J)=0.0 SM000820 

WRITE (6,3) SM000830 

WRITE (6,2) ( (DIFF(I, J) , Jsl ,16), 1=1 ,NTIMES) SM000840 

WRITE (6,2) ( (DIFF(I,J) ,Js17,NTIMES) ,Is1 ,NTIMES) SM000850 

WRITE (16, 4) 3M000860 

FORMAT (’GRID SMOOTHED, ANALYSIS BEGUN *) SM000870 

MIN s 0 SMOOO 88 O 

INT s 0 SM000890 

SCALE s 0.0 SM000900 

CALL BONTUR (P.ISIZE, JSIZE, MIN, INT, SCALE) SM000910 

MIN s 0 SM000920 

INT = 1 SM000930 

SCALE s 1.0 SM000940 

CALL BONTUR (U,ISIZE, JSIZE, MIN, INT, SCALE) SM000S30 

SM000960 

=== BEGINNING OF READ HIGH DENSITY MAP DATA ==========================SM000970 

SM000980 

DO 609 JJJsl, JSIZE SM000990 

DO 609 111=1 , ISIZE SM001000 

609 P(III,JJJ)=-1.0 SM001010 

CALL IN(XLAT, XLONG, TMAX, TMIN, PREC, L, MAXLAT, MINLAT, SM001020 

•MAXLNG, MINLNG) SM001030 
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SMOOIOUO 
SM001050 
SM001060 
SMOO 1070 
SMOO 1080 
SMOO 1090 
SM001100 
sMOomo 
SHOO 1120 
SM001130 
SMOOI 140 
SM001150 
SMOOI 160 


==== BEGINNING OP INTEGRATION OVER A POLYGON ON THE MAP ====s===s=s===sSM001 170 

SMOOI 180 

READ (5,10B,END=9200) IPOLY SM001190 

DO 268 INUM=1 f IPOLY SMOOI 200 

IBUF(2)=INUM SM001210 

READ (5,108,END=9200) L,IDEBUG SMOOI 220 

108 FORMAT (212) SMOOI 230 

N0s1+L*2 SM091240 

M=3 SMOOI 250 

702 READ (5,100,END=9200) XLAT(M) ,XLONG(M) SM001260 

100 FORMAT (F9.3, 1 »F10.3,A5) SMO01270 

XI = ( MAXLAT-XLAT (M) ) *ISCAIE+ 1 SMOO 1 280 

XJ= ( MAXLNG-XLONG (M ) ) • JSCALE+ 1 SMOO 1 290 

I=XI SMOO 1300 

J=XJ SMOOI 310 

IF (XI-FLOAT(I) .GE.0.5) 1=1+1 SM001320 

IF (XJ-FLOAT(J). GE.0.5) J=J+1 SMOOI 330 

IBUF(M)=J SMOO 1340 

IBUF(M+1 )=I SMOOI 350 

IF (M.EQ.3) GO TO 703 SM001360 

IF (IBUF(M) .NE.IBUF(M-2) .0R.IBUF(M+1 ).NE.IBUF(M-1 ) ) GO TO 703 SM001370 

M=M-2 SMOOI 380 

N0=N0-2 SMOOI 390 

703 IF (M.GE.NO) GO TO 701 SM001400 

M=M+2 SMOOI 410 

GO TO 702 SMOOI 420 

701 CALL POLYG (NO , IBUF , IDEBUG f 1X6 ,J , IYMIN , IYMAX ) SMOOI 430 

SMOO 1440 

==== PAINTING OF INTERIOR OF POLYGON WILL COMMENCE ====================SM001450 

SMOOI 460 

I=IYMIN SMOOI 470 

SUM=0.0 SMOOI 480 

ICNT=0 SMOOI 490 

SSsO.O SMOO 1500 

SSIsO.O SM001510 

IC0UNT=0 SMOO 1520 

K=1 SMOOI 530 

206 CONTINUE SMOO 1540 

K1=IX6(K) SMOO 1550 

IF (IX6(K) .GT.5000) K1=K1-5000 SM001560 

K2=IX6(K+1 ) SMOO 1570 

IF(K2.GT.JSIZE) WRITE (16,2RS) SM001580 


IF (L.LT.5) GO TO 503 
DO 502 M=1 ,L 

XI s (MAXLAT-XLAT (M ) ) *ISCALE+ 1 
XJ=(MAXLNG-XLONG(M) )*JSCALE+1 
I=XI 

IF (XI-FLOAT(I) .GE.0.50) 1=1+1 
J=XJ 

IF(XJ-FLOAT(I). GE.0.5) J=J+1 
IF (K.EQ.2) P(I,J)=TMAX(M) 

IF (K.EQ.O) P(I,J)=TMIN(M) 

502 IF (K.EQ.1) P(I,J)=PREC(M) 

503 CONTINUE 
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286 FORMAT (» POLYGON EXTENDS OUTSIDE OBJECTIVE FIELD IN LONGITUDE ') SM001590 
IF (K2.GT.JSIZE) CALL EXIT SM001600 

DO 107 JJ=K1 ,K2 SM001610 

IC0UNT=IC0UNT+1 SM001620 

IF(L.LT.5) GO TO 107 SM001630 

IF(P(I,JJ).LT.-1.0E-10) GO TO 107 SM001640 

ICNT=ICNT+1 SM001650 

SS=SS+U(I,JJ)-P(I,JJ) SMOOI 66 O 

SS1sSS1+(U(I, JJ)-P(I, JJ))**2 SM001670 

107 SUMsSUM+U(I , JJ) SM001680 

IF(I.GT.ISIZE) WRITE (16,296) SM001690 

296 FORMAT ( » POLYGON EXTENDS OUTSIDE OBJECTIVE FIELD IN LATITUDE ' ) SM001700 
IF (I.GT.ISIZE) CALL EXIT SM001710 

IF (IXf>(K+2) .GT.5000) 1=1+1 SM001720 

K=K+2 SM001730 

IF(K.LT.J) GO TO 206 SM001740 

SUM=SUM/ FLOAT ( ICOUNT ) SM001750 

WRITE(7,267) IBUF(2) ,SUM,ICNT SM001760- 

267 FORMAT (• AVERAGE OVER AREA »,I4,’ EQUALS *,F9.5,2X, SM001770 

•’NUMBER OF 0BS=',I5) SM001780 

IF(L.LT. 5) GO TO 268 SM001790 

SM=SS/ ( FLOAT (ICNT)) SM001800 

WRITE (7,362) SM SM001810 

362 FORMAT (* SAMPLE BIAS r F9.5) SM001820 

SMsDSQRT ( ( SS 1 - ( SS« *2 ) / ( DFLOAT ( ICNT ) ) ) / ( DFLOAT ( ICNT- 1 ) ) ) SM001 830 

SSI =SS1/( FLOAT ( ICi.T) ) SM001840 

WRITE (7,363) SM,SS1 SM001850 

363 FORMAT (» STANDARD DEVIATION r *, F9.5, ’ MSE s * F10.5,//) SM001860 

268 CONTINUE SM001870 

IF (L.LT.5) GO TO 5 SM001880 

GO TO S SM001890 

9200 STOP SM001900 

AiND SM001910 

SUBROUTINE IN ( XLAT , XLONG , TMAX , TMIN , PREC , L , MAXLAT , MINL AT , SM001920 

•MAXLNG , MINLNG ) SM001930 

IMPLICIT INTEGER*4 (A-Q, S-Z), REAL* 8 (R) SM001940 

REAL«8 XLAT (500) ,XLONG(500) , TMAX (500) , TMIN ( 500 ), PREC (500) SM001950 

REAL*8 MAXLAT , MINL AT , MAXLNG , MINLNG SM001960 

L=0 SM001970 

' L=L+1 SM001980 

READ(5,100,END=9200)XLAT(L),XLONG(L),TMAX(L),TMIN(L),PREC(L),AAA, SM001990 
* AAB SM002000 

100 FORMAT (F8.3,lJF10.3,2A4) SM002010 

IF (XLAT(L) .LT.-90.0) GO TO 40 SM002020 

ITESTrO SM002030 

IF ( XL AT ( L ) . GT . MAXLAT . OR . XLAT (L).LT. MIWLAT . OR . XLONG (L ) . GT . MAXLNG SM002040 
1 . OR . XLONG ( L ) . LT .MINLNG ) ITESTsI SM002050 

IF (ITEST.EQ.1 ) L=L-1 SM002060 

GO TO 1 SM002070 

40 CONTINUE SM002080 

L=L-1 SM002090 

WRITE (7,106) L SM002100 

’06 FORMAT (• NUMBER OF STATIONS READ IN s *,110,//) SM002110 

9200 RETURN SM002120 

END SM002130 
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POLYG FORTRAN 

WRITTEN BY DAVE PITTS (J.5.C.) 

MODIFIED BY KEVIN MCCULLEN (L.A.R.S.) 06/09/80 

POLYG IS A SUBROUTINE TO THE MAINV WEATHER DATA SMOOTHING 
ALGORITHIM 


VARIABLES 

NO 

BUF 

Xu 

J 

BUF( 1 ) 
BUF(2) 


LENGTH OF INPUT, BUF 
INPUT DATA STRING 
OUTPUT STRING 
LENGTH OF OUPUT, X6 
UNUSED 

NUMBER IDENTIFYING CRD OR PSEUDOZONE 


NO s 1 + L * 2 (L = NUMBER OF PAIRS OF DATA POINTS) 
IF IDEBUG IS 1 THEN EXTRA DEBUG OUTPUT IS GENERATED 


SUBROUTINE POLYG (INO,IBUF , IDEBUG ,1X6 ,IJ ,IYMIN ,IYMAX) 

IMPLICIT INTEGER (A-Q, S-Z), REAL*8 (R) 

INTEGER* H IBUF ( 80 ) , 1X6 ( 5 1 2 ) 

COMMON /STUFF/X6 ( 5 1 2 ) , NPRT , BUF ( 80 ) ,N0,X1 (50) ,Y1 (50) ,N1 ,YMIN,YMAX, 
•X2(55),Y2(55),N2,X3(70),Y3(70),N3,XJ|(512),YlK512),Ni<,X5(200,11),J 
NPRT=16 
NOsINO 
DO 30 1=1 ,80 

30 BUF(I)=IBUF(I) 

102 FORMAT (214) 

CALL SOI (IDEBUG) 

CALL SI 2 (IDEBUG) 

CALL S23 (IDEBUG) 

CALL S3 1 * (IDEBUG) 

CALL SlJ5( IDEBUG) 

CALL S55( IDEBUG) 

CALL S56( IDEBUG) 

IJsJ 

DO 31 1=1,512 

31 IX6(I)=X6(I) 

IYMINsYMIN 

IYMAXsYMAX 

RETURN 

END 


SMS00010 
:SMS00020 
SMS00030 
SMS00040 
SMS00050 
SMS00060 
SMS00070 
SMS00080 
SMS00090 
SMS00100 
:SMS001 10 
SMS00120 
SMS00130 
SMS00140 
SMS00150 
SMS00160 
SMS00170 
SMS00180 
SMS00190 
SMS00200 
SMS00210 
SMS00220 
SMS00230 
SMS00240 
:SMS00250 
SMS00260 
SMS00270 
SMS00280 
SMS00290 
SMS00300 
SMS00310 
SMS00320 
SMS00330 
SMS00340 
SMS00350 
SMS00360 
SMS00370 
SMS00380 
SMS00390 
SMS00400 
sMsooino 

SMSQ0420 

SMS0Q430 

SMS00440 

SMS00450 

SMS00460 

SMS00470 

SMS00480 

SMS00490 

SMS00500 


oooooooooooooook^ooooaooooooooooooooooooooo 
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c 5MS00510 

CBSBBBBSSBSSSSSB3SBSSSSBBEBSBEB8SSSSSXSSSSS*Br BSSBBSEBSBSSBSSSBS*B»®BSB*SSMS00520 

5MS00530 

SOI FORTRAN SMS00540 

WRITTEN BY DAVE FITTS (J.S.C.) SMS00550 

MODIFIED BY KEVIN MCCULLEN (L.A.R.S.) 06/09/60 SM200560 

SMS005 10 

SOI IS A SUBROUTINE IN THE MAINV WEATHER SMOOTHING ALGORITHIM SMSOO50O 

SMS00590 

BBSSSBSBBJ!SSSSSSSB'SSSSSBSSSSSSSSSSBBBBBBBSSSBBBSSSSBSBSSSSSaSSBBSSSBgBSMS00600 

SMS00610 

DEBUG SHOULD BE SET TO 1 TO GET EXTRA PRINTOUT OF WORKING ARRAYS SMS00620 
INPUT IS THRU DATA STRING BUF CONTAINING X AND Y COORDINATES SMS00630 
READS BUF INTO XI AND Y1 AND FINDS YMIN AND YMAX SMS00640 

BUF (I) CONTAINS X(I) IN ODD I STARTING WITH 3 SMS00650 

THE FOLLOWING IS AN EXAMPLE RUN WITH INPUT SMS00660 

1 3 SMS00670 

2 1 SMS00680 

4 3 SMS00690 

6 1 SMS00700 

5 6 SMS00710 

WITH OUTPUT AS FOLLOWS (FIRST LINE IS FOR YMIN) SMS00720 

(LAST LINE IS FOR YMAX) SMS00730 

5002 2 6 6 SMS00740 

5002 3 5 6 SMS00750 

5001 1| 14 6 SMS00760 

5002 5 SMS00770 

5004 5 SMS00780 

5005 5 SMS00790 

BUF(I) CONTAINS Y(I) IN EVEN I STARTING WITH 4, X(3),Y(4) IS A SMS00800 

PAIR, BUF(2) IS A FIELD OR PSEUDOZONE NUMBER, BUF(1) IS BLANK. SMS00810 

OUTPUT IS A DATA STRING X 6 FROM SUBROUTINE S56 SMS00820 

TYPICAL CALLING ROUTINE IS AS FOLLOWS SMSO083O 

CALL SOI SMS00840 

CALL SI 2 SMS00850 

CALL S23 SMS00860 

CALL S34 SMS00870 

CALL S45 SMS00880 

CALL S55 SMS00890 

CALL S56 SMS00900 

BUF AND X 6 ARE IMPLICIT INTEGERS SMS00910 

SMS00920 

SMS00940 

SUBROUTINE SOI (DEBUG) SMS00950 

IMPLICIT INTEGERS (A-Q, S-Z), REAL *8 (R) SMS00960 

COMMON /STUFF/X 6 (512), NPRT , BUF ( 80 ) , NO , X 1 ( 50 ) , Y1 ( 50 ) , N 1 , YMIN , YMAX , SMS00970 
•X2(55),Y2(55),N2,X3(70),Y3(70),N3,X4(512),Yl*(512),N4,X5(200 f 11),J SMS00980 
YMIN= 10000 SMS00990 

YMAXr-IOOOS SMS01000 

j -1 SMS01010 


o o o 
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DO 1 1*3, NO, 2 3M501U2U 

XkBUF(I) SMS01030 

y*BUF(I^1) SMSOIO^Q 

IF (Y.LT. YMIN) YMINsY SMS01050 

IF (Y.GT.YMAX) YMAXsY SMS01060 

J»J+1 SMS01070 

X1(J)*X SMS01080 

1 Y1(J)sY SMSO109O 

NIsJ SMSOIIOO 

SMS01110 


sss* MAKE OUTLINE OVERLAP AT BOTH ENDS sssssssssssssssmsssssbsssssssssSMSOI 120 


102 


XI(l)sXKJ) 

Y1(1)sY1(J) 

J*J+1 

XI (J)sX1 (2) 

Y1(J)=Y1(2) 

IF (J.CT.50) WRITE (16,102) BUF(2) 
FORMAT (* FIELD *,I5,* EXCEEDS THE 
RETURN 
END 


SIZE ALLOWED FOR XI AND 


SMS01130 

SMSOIIOO 

SMS01150 

SMS01160 

SMS01170 

SMS01180 

SMS01190 

Y1»)SMS01200 

SMS01210 

SMS01220 


UUOOODUOUOOOO udu 060 uuo 


3fl 


SSSSSggsSSKKKKXSSSgSSSSSSSS; 


SB'fESS 


SI 2 FORTRAN 

WRITTEN BY DAVE FITTS (J.S.C.) 
MODIFIED BY KEVIN MCCULLEN (L.A.R.S.) 


06/09/80 


S12 FORTRAN IS PART OF MAINV WEATHER DATA SMOOTHING ALGORITHIM, 
SI 2 REMOVES REDUNDANT POINTS FROM XI AND Y1 SO THAT THERE ARE 
AT MOST 2 CONTIGUOUS POINTS ON A LINE. 


:ss:sssr 


SMS 01230 

:ce7«ssssessesssxbssssssSMS01240 

SMS01250 
SMS01260 
SMS 01270 
SMS012B0 
SMSO1290 
SMS01300 
SMS01310 
SMS01320 
SMS01330 
:SMS013*!0 
SMS01350 
SMS01360 
SMS01370 
SMS013B0 


SUBROUTINE SI 2 (DEBUG) 

IMPLICIT INTEGER**! (A-Q f S-Z), REAL*8 (R) 

COMMON /STUFF/X6(i>12) ,NPRT,BUF(80) ,N0,X1 (50) ,Y1 (50) ,N1 ,YMIN,YMAX, 
•X2(55),Y2(55),N2 f X3(70),Y3(70) l N3 f Xl!(512),Y*!(512),Nl! > X5(200,11),J SMS01390 
J=1 SMS01*!00 

DO 1 ls2,N1 SMS01410 

IF (YI(I).NE.YI(I-l)) GO TO 2 SMS01*!20 

IF (Y1(I).NE.Y1(I+1)) GO TO 2 SMS01H30 

SMS01440 

POINT I IS A REDUNDANT POINT =r========= = r- = ================B=== = sESMS01 *!50 

SMS0 1*160 
SMS0 1*170 
SMS01*JB0 
SMS01490 

:=== POINT IS NOT A REDUNDANT POINT r========= = === = ==== = == == = == ===== = ==SMS01500 

SMS01510 
SMS01520 
SMS01530 
SMS01540 

1 CONTINUE SMS01550 

SMS01560 
SMS01570 

===== MAKE OUTLINE OVERLAP AT BOTH ENDS ====================b===========SMS01580 

SMS01590 

X2(1 )=X2(J) SMS01600 

Y2(1 )sY2(J) SMS01610 

J=J*1 SMS01620 

X2(J)sX2(2) SMS01630 

Y2( J)=Y2(2) SMS01640 

JsJ-fl SMS01650 

X2(J)=X2(3) SMS01660 

Y2(J)*Y2(3) SMS01670 

IF(J.GT.55) WRITE (NPRT,102) BUF(2) SMS01680 

102 FORMAT (» FIELD * ,15,* EXCEEDS THE SIZE ALLOWED FOR X2 AND Y2')SMS01690 
RETURN 5MS01700 

END SMS01710 


GO TO 1 
CONTINUE 


J=J+1 

X2(J)=X1 (I) 
Y2( J)=Y1 (I) 
CONTINUE 
N2=J 
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C SMS01720 

Csszxsscsssssbss=sss==sss=szsss=ssssbssssbs=sssssssesss==sss5:ss:=:s=ss:SMS01730 
C SMS01740 

C S23 FORTRAN SMS01750 

C WRITTEN BY DAVE PITTS (J.S.C.) SMS01760 

C MODIFIED BY KEVIN MCCULLEN (L.A.R.S.) 06/09/80 SMS01770 

C SMS01780 

C S23 IS PART OF THE MAINV WEATHER DATA SMOOTHING ALGORITHIM, S23 SMS01790 

C INSERTS REDUNDANT POINTS AT MAXIMA, MINIMA, AND INFLECTION POINTS SMS01800 

C SMS01810 

C== = = === = = = = = = = = = = = - = = = = = = = = = = = = = = r::=== = = = = === = = = = = = = = === = = = = = .- = = = = = = ==SMS0l820 
C SMS01830 

SUBROUTINE S23(DEBUG) SMS01840 

IMPLICIT INTEGER* 1 ! (A-Q, S-Z), REAL*8 (R) SMS01850 

COMMON /STUFF/X6(512) ,NPRT,BUF(80) , NO, XI (50) ,Y1 (50) ,N1 ,YMIN,YMAX, SMS01860 
•X2(55) ,Y2(55) ,N2,X3(70) ,Y3(70) ,N3 ,Xi| (512) , Y4(512) ,N4,X5(200,1 1 ) , J SMS01870 
JzO SMfi 31880 

DO 1 Is2,N2 SMS01890 

YY1=Y2(I-1 ) SMS01900 

YY2=Y2(I ) SMS01910 

YY3*iY2(I+1 ) SMS01920 

YY4bY2(I+2) SMS01930 

D12=YY2-YY1 SMS01940 

D23=YY3-YY2 SMS01950 

C SMS01960 

C==== CHECK TO SEE IF POINTS I AND (1+1) ARE POINTS OF INFLECTION ==r===SMS01970 

C SMS01980 

IF(D23.EQ.O) GO TO 2 SMS01990 

C SMS02000 

C==== CHECK TO SEE IF POINTS I ArtD (1-1) ARE A TWO-POINT MAX OR MIN ====SMS02010 

C SMS02020 

IF(D12.EQ.O) GO TO 3 SMS02030 

C SMS02040 

C==== CHECK TO SEE IF POINTS I AND (1-1) ARE A ONE-POINT MAX OH MIN ====SMS02050 

C SMS02060 

IF ( (D12.GT.0) .AND. (D23.GT.0)) GO TO 3 SMS02070 

IF((D12.LT.O). AND. (D23.LT. 0)) GO TO 3 SMS02080 

C SMS02090 

C --:=== POINT I AIS A MAXIMUM OR MINIMUM =================================SMS02100 

C SMS02110 

J=J+1 SMS02120 

X3(J)=X2(I) SMS02130 

Y3(J)*Y2(I) SMS02140 

GO TO 3 SMS02150 

2 CONTINUE SMS02160 

C SMS02170 

C==== POINTS I AND (1+1) MIGHT BE POINTS OF INFLECTION =================SMS02l80 

C SMS02190 

D34=YY4-YY3 SMS02200 

IF((D12.GT.O).AND.(D34.LT.O)) GO TO 3 SMS02210 

IF( (D12.LT.0) .AND. (D34.GT.O)) GO TO 3 SMS02220 


ooo ooo n n o 


:::« POINTS I AND (1+1) ARE POINTS OF INFLECTION =*====="============= 

J=J+1 

Y3(J)=Y2(I) 

IF(X2(I+1 ).LT.X2(I)) CO TO 4 

==== PUT A REDUNDANT POINT TO RIGHT OF POINT I AND TAG BY ADDING 5000 = 

X3(J)=X2(I)+5001 
GO TO 3 
J* CONTINUE 

»sr PUT A REDUNDANT POINT TO LEFT OF POINT I AND TAG BY ADDING 5000 == 


X3(J)=X2(I)+4j?99 
3 CONTINUE 
JsJ+1 

X3(J)=X2(I) 

Y3(J)=Y2(I) 

1 CONTINUE 
J=J+1 

X3(J)=X3(1 ) 

Y3v'J)=Y3(1) 

N3-J 

IF (J.GT.70) WRITE (NPRT,102) BUF(2) 

102 FORMAT (' FIELD ',15,' EXCEEDS THE SIZE ALLOWED FOR X3 AND Y3V 
RETURN 
e:?d 



o n oooooooocs oooooooooooo 
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S3 1 ! FORTRAN 

WRITTEN BY DAVE PITTS (J.S.C.) 
MODIFIED BY KEVIN MCCULLEN (L.A.R.S.) 


06/09/80 


S3 1 * IS PART OF THE MAINV WEATHER DATA SMOOTHING ALGORITHIM, S3 1 * 
FILLS IN MISSING LINES IN THE DATA SET 


SUBROUTINE S3*J (DEBUG) 

IMPLICIT INTEGERS (A-Q, S-Z), REAL*8 (R) 

COMMON /STUFF/X6(512) ,NPRT,BUF(80) ,N0,X1 (50) ,Y1 (50) ,N1 t YMIN,YMAX f 
•X2(55) ,Y2(55) ,N2,X3(70) ,Y3(70) ,N3,Xl<(512) ,Y^I(512) ,N»I,X5(200, 1 1 ) , J 
N3=N3-1 

J=0 

DO 1 Isl,N3 
XB=X3(I) 

YB=Y3(I) 

XN=X3(I+1) 

YNsY3(I+1) 

J=J+1 

X4(J)=XB 

Y^»(J)sYB 

INC=YN-YB 

IF(INC.EQ.O) GO TO 1 


MISSING LINES MUST BE FILLED IN. 

CHECK TO SEE IF EITHER I OR (1+1) HAS BEEN TAGGED AS POINT OF 
INFLECTION. 


IF(XB.GT.3000) XB=XB-5000 
IF(XN.GT.3000) XN=XN-5000 
RDXsDFLOAT(XN-XB) 

RDY =DFLOAT ( YN -YB ) 

RS=RDX/RDY 

RXB=DFLOAT(XB)+0.5 

INC=1 

IF(RDY.LT.O.O) INC=-1 
Y--YB 

FILL IN LINES BETWEEN (BUT NOT INCLUDING) POINTS I AND (1+1) 


3 CONTINUE 


SKS02520 

=SMS02530 

SMS02540 

SMS02550 

SMS02560 

SMS02570 

SMS02580 

SMS02590 

SMS02600 

SMS02610 

=SMS02620 

SMS02630 

3MS0261J0 

SMS02650 

SMS02660 

SMS02670 

SMS02680 

SMS02690 

SMS02700 

SMS02710 

SMS02720 

SMS02730 

SMS027 1 »0 

SMS02750 

SMS02760 

SMS02770 

SMS02780 

SMS02790 

SMS02800 

SMS02810 

SMS02820 

SMS02830 

SMS02840 

SMS02850 

SMS02860 

:SMS02870 

SMS02880 

SMS02890 

SMS02900 

SMS02910 

SMS02920 

SMS02930 

sMsopg^o 

SMS 02950 

SMS 02960 

SMS 02970 

:SMS02980 

SMS02990 

SMS03000 
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JaJ+1 

YsY+INC 

RYaDFLOAT(Y-YB) 

RXaRXB+RS«RY 

X4(J)=RX 

Y4(J)=Y 

IF(Y.NE.YN) GO TO 3 
JsJ-1 

1 CONTINUE 
NtyaJ 

IF(J.GT.511) WRITE (NPRT.102) BUF(2) 


SMS03010 

SMS03020 

SNS03030 

SMS03040 

SMS03050 

SMS03060 

SMS03070 

SMS03080 

SMS03090 

SMS03100 

SMS03110 


102 FORMAT (• FIELD ',15,' EXCEEDS THE SIZE ALLOWED FOR X4 AND Y4') SMS03120 


RETURN 


SMS03130 


END 


SMS031M0 



nooooooooonrj 
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SMS03150 

ftssassssssssssssBBssssssBSSsssrssssssssrsBssrssBSSBsssssssrsssssr ssftsssSMS03l60 

SMS03170 

S45 FORTRAN 

WRITTEN BY DAVE PITTS (J.S.C.) 

MODIFIED BY KEVIN MCCULLEN (L.A.R.S.) 06/09/60 


S34 IS PART OF THE MAIV WEATHER DATA SMOOTHING ALGORITHIM, S45 
COLLECTS ALL OF THE INTERCEPTS WITHIN GIVEN LINES 


SMS03180 
SMS03190 
SMS03200 
SMS03210 
SMS03220 
SMS03230 
SMS032H0 

BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBr,B=SBBBBBB=BBBBBBBBBBBBBBB=BBB=BBBBBBBBSMS03250 

SMS03260 

SUBROUTINE S45(DEBUG) SMS03270 

IMPLICIT INTEGERS (A-Q, S-Z), REAL* 8 (R) SMS03280 

COMMON /STUFF/X6(512) ,NPRT,BUF(8C) ,N0, XI (50) ,71 (50) ,N1 ,YMIN,YMAX, SMS03290 
•X2(55) ,Y2(55) ,N2,X3(70) ,Y3( 70) ,N3,X4(512) ,Yi|(512) ,N4 ,X5(200, 1 1 ) , J SMS03300 
YOFFsYMIN-1 SMS03310 

YENDrYMAX-YOFF SMS03320 

IF (YEND.GT.P0O) WRITE (NPRT,102) BUF(2)| SMS03330 

102 FORMAT (• FIELD ',15,' HAS TOO MANY LINES' ) SMS03340 

IF (YEND.GT.200) STOP SMS03350 

DO 1 Ib 1,200 SMS03360 

X5(I,11)=0 SMS03370 

1 CONTINUE SMS03380 

IF (iM.GT.512) WRITE (NPRT,200) N4 SMS03390 

200 FORMAT (• Nl» s ',15) SMS03400 

IF (N4.GT.512) STOP SMS03410 

SMS03420 
SMS03430 
SMS03440 
SMS03450 
SMS03460 


DO 2 1=1, N4 
SeX4(I) 

L=Y4(I)-YOFF 

IF (L.GT.200) WRITE (NPRT,201) L 
201 FORMAT (• Lb \I5) 

IF (L.GT.200) STOP 
N=X5(L, 1 1 ) 

N=N+1 

IF (N.GT.10) WRITE (NPRT,103) L,BUF(2) 
103 FORMAT (• -LINE ',15,' OF FIELD »,I5,' 
•NS • ) 

IF (N.GT.10) STOP 
X5(L, 1 1) b N 
X5(L,N)=S 
2 CONTINUE 
DO 3 L=1,YEND 
NERDsX5(L,ri) 


SMS03470 

SMS03480 

SMS03490 

SMS03500 

HAS TOO MANY INTERSECTIOSMS035 1 0 

SMS03520 

SMS03530 

SMS03540 

SMS03550 

SMS03560 

SMS03570 

SMS03580 


3 CONTINUE 
RETURN 
END 


SMS03590 

SMS03600 

SHSQ3610 


ooo o o o onoooooo ooo oooooooooooo 


0-<m 
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SMS03620 
SMS03640 

S55 FORTRAN SMS03650 

WRITTEN BY DAVE PITTS (J.S.C.) SMS03660 

MODIFIED BY KEVIN MCCULLEN (L.A.R.S.) 06/09/80 SMS03670 

SMS03680 

S55 IS PART OF THE MAINV WEATHER DATA SMOOTHING ALOORITHIM, S55 SMS03690 

SORTS THE INTERCEPTS INTO ASCENDING ORDER SMS03700 

SMS03710 

:=====zz================z====z=zzzzz===zz=z=z=====z==zzz=z=====zz==zzz=SMS03720 

SMS03730 

SUBROUTINE S55( DEBUG) SMS03740 

IMPLICIT INTEGER*^ (A-Q, S-Z), REAL* 8 (R) SMS03750 

COMMON /STUFF /X6 ( 5 1 2 ) , NPRT , BUF ( 80 ) ,N0,X1 (50) ,Y1 (50) ,N1 ,YMIN,YMAX, SMS03760 
•X2(55) ,Y2(55) ,N2,X3(70) ,Y3(70) ,N3,X4(512) ,Y4(512) ,N4,X5(200, 1 1 ) , J SMS03770 
YENDs YMAX- YMIN-f 1 SMS03780 

DO 1 L=1,YEND SMS03790 

NEND=X5(L,11) SMS03800 

NODDsNEND-2* ( NEND/2 ) SMSO 38 IO 

IF (NODD.EQ.O) GO TO 6 SMS03820 

SMS03830 

=== AN ODD NUMBER OF INTERSECTIONS IS NOT PERMITTED ==r========= == == == SMS03840 

SMS03850 

LINEsL+YMIN-1 SMS03860 

102 FORMAT ( » ODD NUMBER OF VERTICES ON LINE *,I5,’ OF FIELD • ,15) SMS03870 
6 CONTINUE SMS03880 

DO 2 1=1 ,NEND SMS03890 

XMINs 30000 SMS03900 

DO 3 Jsl ,NEND SMS03910 

X=X5(L,J) SMS03920 

SMS03930 

SMS03950 

IF THE POINT HAS BEEN USED BEFORE (AND TAGGED AS 31000) JUMP OVER SMS03960 
IT SMS03970 

SMS03980 

SMS04000 

IF (X.EQ.31000) GO TO 3 SMS04010 

SMS04020 

=== IF THE POINT IS TAGGED AS POINT OF INFLECTION SUBTACT 5000 =======SMS04030 

SMS04040 

IF (X.GT.3000) X=X-5000 SMS04050 

IF (X.GT.XMIN) GO TO 3 SMS04060 

XMIN=X SMS04070 

JMIN=J SMS04080 

3 CONTINUE SMS04090 

SMS04100 

=== POINT STORED AT JMIN HAS THE SMALLEST REMAINING X-VALUE ==========SMS04l 10 

SMS04120 

XI (I)=X5(L, JMIN) SMS04130 


o o o o o o non 
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C 

C=r=r TAG POINT AT JMIN AS HAVING BEEN USED ================= 

C 

X5 (L , JMIN ) =31000 
2 CONTINUE 
DO 4 1=1 ,NEND 
XsXI(I) 

IF (X.LT.3000) GO TO 5 

==« THIS POINT IS A NECESSARY REDUNDANT POINT OF INFLECTION 

JsI/2 
SW=I-2*J 


===* POINT IN EVEN POSITION SHOULD BE MOVED TO RIGHT ======= 

IF(I.EQ.NEND) GO TO 969 
IF (SW.EQ.O) X5(L,I)=X1(I+1 ) 

==== POINT IN ODD POSITION SHOULD BE MOVED TO LEFT ========= 

IF(I.EQ.I) GO TO 5 

969 IF (SW.NE.O) X5(L,I)=X1(I-1) 

GO TO 4 

5 CONTINUE 
X5(L,I)=X1(I) 

4 CONTINUE 
DO 970 1=1 ,NEND 

IF (X5(L,I) .EQ.31000) X5(L,I)=X1(I) 

IF (X5(L,I).GT.5000) X5(L,I)=X5(L,I)-5000 

970 CONTINUE 
1 CONTINUE 

RETURN 

END 


SMS04140 
:SMSO4l50 
SMS04160 
SMS 04170 
SMS04180 
SMS04190 
SMS04200 
SMS04210 
SMS 04 220 
;SMS04230 
SMS04240 
SMS04250 
SMS04260 
SMS04270 
:SMS04280 
SMS04290 
SMS04300 
SMS04310 
SMS04320 
:SMS04330 
SMS04340 
SMS04350 
SMS04360 
SMS04370 
SMS04380 
SMS04390 
SMS04400 
SMS04410 
SMS04420 
SMS04430 
SMS04440 
SMS04450 
SMS04460 
SMS04470 


uuouuuuouooa ouo 
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::::::::::::::::: 




S56 FORTRAN 

WRITTEN BY DAVE PITTS (J.S.C.) 
MODIFIED ?i KEVIN MCCULLEN (L.A.R.S.) 


06/09/80 


S56 IS PART OF THE MAINV WEATHER DATA SMOOTHING ALGORITHIM, S56 
PACKS THE INTERCEPTS INTO A ONE-DIMENSIONAL BUFFER 


SUBROUTINE S56( DEBUG) 

IMPLICIT INTEGERS (A-Q, S-Z), REAL* 8 (R) 

COMMON /STUFF/X6(512),NPRT,BUF(80) f N0,X1(50),Y1(50),N1,YMIN f YMAX, 
*X2(55) , Y2(55) , N2,X3(70) ,Y3<70) ,N3,X4<512) ,Y4(512) ,N4,X5(200, 1D.J 
YEND=YMAX-YMIN+1 
J=0 

DO 1 L=1,YEND 

=S== tag THE FIRST X-VALUE IN EACH LINE BY ADDING 5000 ================ 

X5(L,1)=X5(L,1)+5000 
NEND=X5(L, 1 1 ) 

905 IQENDrNEND-3 

IF (NEND.LT.4) GO TO 903 
DO 902 I*1,IQEND,2 

IF (X5(L,I+1 ).NE.X5(L,I+2)) GO TO 902 
X5(L,I+1 )=X5(L,I+3) 

IBEG=I+2 
IEND=IQEND+1 
NENDsNEND-2 

IF (IBEG.GT.IQEND) GO TO 902 
DO 904 IX=IBEG f IEND 
904 X5(L,IX)sX5(L,IX+2) 

GO TO 905 

902 CONTINUE 

903 CONTINUE 
DO 2 1=1 ,NEND 
J=J+1 

X6(J)=X5(L,I) 

2 CONTINUE 
1 CONTINUE 

IF (J.GT.511) WRITE (NPRT.102) BUF(2) 

102 FORMAT (• FIELD' ,2X,I5, ' EXCEEDS THE SIZE ALLOWED FOR X6 ') 
N6=J 

X6(512)=N6 
RETURN 
END 


SMS04480 

SM504490 

5MS04500 

SMS04510 

SMS04520 

SMS04530 

SMS04540 

SMS04550 

SMS04560 

SMS04570 

SMS04580 

SMSO4590 

SMS04600 

SMS04610 

SMS04620 

SMS04630 

SMS04640 

SMS04650 

SMS04660 

SMS04670 

SMS04680 

SMS04690 

SMS04700 

SMS04710 

SMS04720 

SMS04730 

SMS04740 

SMS04750 

SMS04760 

SMS04770 

SMS04780 

SMS04790 

SMS04800 

SMS04810 

SMS04820 

SMS04830 

SMS04840 

SMS04850 

SMS04860 

SMS04870 

SMS04680 

SMS04890 

SMS04900 

FMS04910 

SMS04920 

SMS04930 

SMS04940 

SMS04950 

SMS04960 
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c 


SMSO4970 

Css 

:pSSSSSSS5SSSSSSSSSSWSS5SSSS&SSSSSS££SSSSrsSSSSSS5?SS£SSSSSSSS5?££S 

sssSMS04980 

c 


SMSO4990 

c 

BONTUR FORTRAN 

SMS05000 

c 

WRITTEN BY DAVE PITTS (J.S.C.) 

SMS05010 

c 

MODIFIED BY KEVIN MCCULLEN (L.A.R.S.) 06/09/80 

SMS05020 

c 


SMS05030 

c 

BONTUR IS PART OF THE MAINV WEATHER DATA SMOOTHING ALGORITHIM, 

SMS05040 

c 

BONTUR PRINTS A STANDARD NI BY NJ GRID WITH CONTURING BETWEEN 

SMS05050 

c 

LINES 

SMS05060 

c 


SMS05070 

Css 

SSSSSSSSSSSSSSSSSSSSZSZSSSSSSZZSSSSSSSrsSSSSSZSSZSSSSSSSSSSSSSSSSS 

ss=SMS05080 

c 


SMS05090 

c 

VARIABLES 

SMS05100 

c 


SMS051 10 

c 

MIN MINIMUM VALUE 

SMS05120 

c 

INT CONTOURING INTERVAL 

SMS05130 

c 

SCALE SCALING FACTOR FOR PRINTING 

SMS05140 

c 


SMS05150 

c 

IF INT s 0 THEN THERE WILL BE NO CONTOURS OR DATA PRINTED 

SMS05160 

c 

IF NJJ IS GREATER THAN 26 2 GRIDS ARE PRINTED; 

SMS05170 

c 

FROM 1 TO 26 

SMS05180 

c 

FROM 26 TO NJJ 

SMS05190 

c 


SMS05200 

Css 


sssSMS05210 

c 


SMS05220 


SUBROUTINE BONTUR (Z,NI,NJJ,MIN,INT, SCALE ) 

SMS05230 


IMPLICIT INTEGERS (I-N), REAL»8(A-H, 0-Z) 

SMS05240 


INTEGERS IZ(64,64) 

SMS05250 


INTEGERS KALP(16) ,LINE(127) ,LIN(27) 

SMS05260 


REAL*8 Z(64,64) 

SMS05270 


DATA KALP/1H ,1HA,1H ,1HB,1H ,1HC,1H ,1HD,1H ,1HE,1H ,1HF,1H , 

SMS05280 


1 1HG.1H , 1HH/ 

SMS05290 


LT0T=INT*16 

SMS05300 


NTEMP s NJJ 

SMS05310 


NJJ s 51 

SMS05320 


NJsNJJ 

SMS05330 

c 


SMS05340 

Css 

ss 360 ss======s=====z===ss===ss====s==s=ss=ss=s===s=rss=s===s==== 

ss=SMS05350 

c 


SMS05360 


Jlsl 

SMSO5370 


IF(NJJ.GT.26) NJs26 

SMS05380 


DO 10 Isl.NI 

SMS05390 


DO 10 Jsl.NJJ 

SMS05400 


10 IZ(I,J)sZ(I,J)»SCALE 

SMS05410 


IF (INT) 51,50,51 

SMS05420 


51 NIMsNI-1 

SMS05430 


60 NJMsNJ-JI 

SMS05440 


WRITE(6,910) 

SMS05450 

910 

FORMAT(IHI) 

SMS05460 


NUM=5*NJM+1 

SMS05470 


WRITE(6,900)(IZ(1,J),JsJ1,NJ) 

SMS05480 

900 

FORMAT (3X, 2615) 

SMS05490 


- •«* *n». 
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DO 1 IRs2,NI 

SMS05500 

DO 2 JDrl t 2 

SMS05510 

IF(JI.NE.I) GO TO 20 

SMS05520 

DO 3 Lrl ,NJ 

SMS05530 

3 LIN(LM<IZ<IR,L)-IZ<IR-1 f L))»JD)/3+IZ(IR-1,L) 

SMS05540 

GO TO 30 

SMS05550 

20 DO 1(0 Ls26,NJJ 

SMS05560 

KO LIN(L-25)=(<IZ<IR,L)-IZ(IR-1 f L))»JD)/3+IZ(IR-1,L) 

SMS05570 

30 K*1 

SMS055B0 

DO 1( J=1,NJM 

SMS05590 

LINJsLIN(J) 

SMS05600 

LINE(K)=LINJ 

SMS05610 

NDZ=LIN(J+1)-LINJ 

SMS05620 

DO 5 L=1 f 4 

SMS05630 

KsK+1 

SMS05640 

5 LINE (K)s(NDZ*L) /5+LIN J 

SMS05650 

K=K+1 

SMS05660 

4 CONTINUE 

SMS05670 

LINE(K)sLIN(NJM+1 ) 

SMS05680 

DO 6 L=1 ,NUM 

SMS05690 

JDFsLINE ( L) -MIN 

SMS05700 

IF(JDF)8,9,9 

SMS05710 

8 JDF=JDF-LT0T* ( ( JDF+1 ) /LTOT- 1 ) 

SMS05720 

9 J=JDF/INT 

SMS05730 

IF(J-16)6, 26,26 

SMS05740 

26 J*J-(J/16)»16 

SMS05750 

6 LINE(L)=KALP(J+1 ) 

SMS05760 

WRITE(6, 901) (LINE(L) ,Ls1 ,NUM) 

SMS05770 

901 FORMAT ( 7X 126A1 ) 

SMS05780 

2 CONTINUE 

SMS05790 

WRITE(6 l 900)(IZ(IR l J) f JsJI ,NJ) 

SKS05800 

1 CONTINUE 

SMS05810 

IF(NJ.NE.NJJ) GO TO 223*1 

SMS05820 

NJJ s NTEMP 

SMS05830 

RETURN 

SMS05840 

2234 CONTINUE 

SMS05850 

NJrNJJ 

SMS05860 

J1=26 

SMS05870 

GO TO 60 

SMS05880 

50 CONTINUE 

SMS05890 

IF (NJ .NE.NJJ) GO TO 2235 

SMS05900 

NJJ = NTEMP 

SMS05910 

RETURN 

SMS05920 

2235 CONTINUE 

SMS05930 

NJsNJJ 

SMS059K0 

J1=26 

SMS05950 

GO TO 50 

SMS05960 

END 

SMS05970 


000000000000000 0 00000 0 0000000000 0 0-0 0000000000000000 
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ANAL FORTRAN 

WRITTEN BY DAVE PITTS (J.S.C.) 

MODIFIED BY KEVIN MCCULLEN (L.A.R.S.) 06/09/80 

ANAL IS PART OF THE MAINV WEATHER DATA SMOOTHING ALGOR ITHIM, ANAL 
PERFORMS VARIATIONAL ANALYSIS BY THE 'KIT WAGNER 2ND DERIVATIVE 
FILTERING* METHOD. 

ANAL MINIMIZES THE INTEGRAL... 

SQUARES OF THE DIFFERENCES + THS SQUARES OF THE GRADIENT* ALF2 4 
THE SQUARES OF THE LAPLACIAN*ALF4 


VARIABLES 

UO INPUT DATA 

U ANALYSIS 

AA FILTER WEIGHTS (SEE NOTE BELOW) 

ALF2 FILTER WEIGHT 

ALF4 FILTER WEIGHT 

MAXPAS MAXIMUM NUMBER OF ITERATIONS 

ERR APROXIMATION CRITERIA 

ARRAY DIMENSIONS 

U(NI,NJ) 

UO(NI,NJ) 

WA(NI,NJ) 

YA(NI+2,NJ+2) 

NOTE: FOR FILTER WEIGHTS REFERENCE K. WAGNER THESIS 


FOR MESOSCALE ANALYSIS OF MAGNITUTE 10, TYPICAL PARAMETERS ARE: 
AA = 100.0 
ALF2 = 1.0 
ALF4 = 1.0 
MAXPAS = 99 
ERR = .01 

INCREASING ALF2 AND/OR ALF4 REDUCES HIGHER FREQUENCYS 

TYPICAL MAXIMUMS ARE: 

ALF2 = 10.0 
ALF4 = 10.0 
AA = 100.0 

VALUES FOR ALF2 AND ALF4 ARE USUALLY .1, 1.0, 10.0 
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SUBROUTINE ANAL ( NI , N J , AA , ALF2 , ALF4 ,MAXPAS , ERR , UO , U ) SMS06U9Q 

IMPLICIT INTEGER - *! (I-N), REAL* 8 <A-H,0-Z) SMS06500 

DIMENSION U(NI,NJ),UO(NI,NJ),YA(66,66),WA(66,66) SMS06510 

EQUIVALENCE (YA( 1 , 1 ) ,WA< 1 , 1 > ) SMS06520 

10*16 SMS0653O 

NJP2sNJ+2 SMS06540 

NIP2xNU2 SMS06550 

NIM1=NI-1 SMS06560 

NIM2=NI-2 SMS06570 

NJMUNJ-1 SMS06580 

NJM2=NJ-2 SMS06590 

NJP1=NJ+1 SMS06600 

NIP1=NI+1 SMS06610 

BETA*2.0 SMS06620 

C SMS06630 

Cxsxx INITIALIZE GUESS FIELD BY AVERAGING =======r=v»===r===*=====^=====sSMS066 i !0 

C SMS06650 

DO 16 Jxl ,NJ SMS06660 

DO 16 Ixl, NI SMS06670 

16 U(I,J)=UO(I,J) SMS06680 

DO 10 Js1,NJP2 SMS06690 

DO 10 1=1 ,NIP2 SMS06700 

10 YA(I,J)=0.0 SMS06710 

DO 9997 J*1,NJ SMS06720 

DO 9997 I*1,NI SMS06730 

IF (U(I,J) .NE. 0) GO TO 9998 SMS06740 

9997 CONTINUE SMS06750 

DO 9996 J=1 ,NJ SMS06760 

DO 9996 1=1 ,NI SMS06770 

U0(I t J)=0.0 SMS06780 

9996 U(I,J)=0.0 SMS06790 

RETURN SMS06800 

9998 KNT=1 SMS06810 

201 CONTINUE SMS06820 

C SMS06830 

C==== CHECK FOR NUMBER OF NO GUESS ==============s=s=xxx=x==sxsx=xxxxxxxSMS06840 

C SMS06850 


IF (KNT) 15,200,15 
15 KNTxO 

DO 12 J=2,NJP1 
DO 12 1=2,NIP1 
12 YA(I,J)=U(I-1 ,J-1 ) 

DO 99 J=2,NJP1 
. DO 99 I=2,NIP1 

IF (YA(I,J) ) 86,98,86 
98 SUMsO.O 
CNTxO.O 


SMS06860 

SMS06870 

SMS06880 

SMS06890 

SMS06900 

SMS06910 

SMS06920 

SMS06930 

SMS06940 

SMS06950 

SMSO 696 O 


Css sx AVERAGE NINE POINTS ======s===s=s=s===s==ssxxxxxxxxxxxxxxxxxrx==xxSMS06970 

C SMS06980 


o u o 
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24 DO 97 JK*1,3 
DO 97 IK* 1,3 
II.I-2+IK 
JJsJ-2+JK 

IF (YA(II,JJ)) 96,97,96 

96 5UM»SUM+YA(II,JJ) 

CNTsCNT+1 . 

97 CONTINUE 

IF (CNT) 93,92,93 

93 IF (SUM) 95,911,95 
C 

Cxssr USE .0001 INSTEAD OF ZERO AVERAGE sssssm8b*s»»»mss»k*sssx: 
C 

94 U(I-l,J-1)r.0001 
GO TO 99 

95 U(I-1 , J— 1 )=SUM/CNT 
GO TO 99 

92 KNTrKNT+1 
GO TO 99 

86 U(I-1 , J-1 )=YA(I, J) 

99 CONTINUE 
200 CONTINUE 

WRITE(IO, 100) KNT 

100 FORMAT ( 15, 29H POINTS UNSPECIFIED THIS PASS) 

IF (KNT) 201,79,201 
79 CONTINUE 

ass SMOOTH FIELD OF AVERAGES -==================================-- 

DO 31 Js2,NJM1 
DO 31 Is2,NIM1 

31 WA(I,J)a(4»*U(I, J)+U(I-1 , J-1 )+U(I+1 ,J-1 )+U(I+1 ,J+1 )+U(I-1 , J+1) 
1+2. # (U(I,J-1 )+U(I+1 ,J)+U(I,J+1)+U(I-1 ,J)))/l6. 

DO 32 Ja2,NJM1 


SMS06990 
SMS07000 
SMS07010 
SMS07020 
SMS07030 
SMS07040 
SMS07050 
SMS07060 
SMS07070 
SMS07080 
SMS07090 
SMS07100 
SMS07110 
SMS07120 
SMS07 1 30 
SMS07140 
SMS07150 
SMS07160 
SMS07170 
SMS07180 
SMS07190 
SMS07200 
SMS07210 
SMS07220 
SMS07230 
SMS07240 
SMS07250 
:SMS07260 
SMS07270 
SMS07280 
SMS07290 
SMS07300 
SMS07310 
SMS07320 


WA(1,J)s(8.*U(1 , J)+2.1(U(2, J)+U(1 ,J-1 )+U(1 ,J+1 ))+U(2, J-1 )+U(2,J+1 SMS07330 
1))/16. SMS07340 

32 WA(NI, J)=(8.*U(NI, J)+2.*(U(NIM1 ,J)+U(NI, J-1 )+U(NI, J+1 ) )+U(NIM1 ,J-1SMS07350 

1 )+U(NIM1 , J+1 ))/l6. SMS07360 

DO 33 Ia2,NIM1 SMS07370 

WA(I, 1 )=(8.*U(I, 1 )+2.*(U(I,2)+U(I-1 , 1 )+U(I+1 ,1 ))+U(I-1 ,2)+U(I+1 ,2 SMS07380 
1))/16, SMS07390 

33 WA(I,NJ)=(8.»U(I,NJ)+2.»(U(I,NJM1)+U(I-1,NJ)+U(I+1,NJ))+U(I-1,NJM1SMS07400 

1)+U(I+1 ,NJM1))/l6. SMS07410 

WA(1,NJ)s(3.*U(1,NJ)+2.»(U(2,NJ)+U(1,NJM1))+U(2,NJM1))/8. SMS07420 

WA(NI,NJ)=(3. # U(NI,NJ)+2.«(U(NIM1,NJ)+U(NI,NJM1))+U(NIM1,NJM1))/8.SMS07430 
WA(NI, 1 )s(3.»U(NI, 1 )+2.*(U(NIM1 , 1 )+U(NI,2) )+U(NIM1 ,2))/8. SMS07440 

WA(l,1)a(3.»U(1,1)+2.*(U(2,1)+U(1,2))+U(2,2))/8. SMS07450 

DO 36 J=1 ,NJ SMS07460 

DO 36 Is1,NI SMSO7470 

36 U(I,J)=WA(I,J) SMS07480 

39 CONTINUE SMS07490 
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c 

SMS07500 


Cxxss 2 DIMENSIONAL ANALYSIS OF INTERIOR POINTS AND WEIGHTS 

sssbbbSMS07510 


C 

SMS07520 


GRDs(NI-4)»(N<7-4) 

SMS07530 


alf4Bxalf4»beta 

SMS07540 


UIJ0x4.«ALF2+20.»ALF4 

SMS07550 


UI1J1b-ALF2-8.»ALF4 

SMS07560 


WRITE (7, 5 00) AA,ALF2,ALF4,UIJ0,UI1J1 

SMS07570 

a 

C 

SMS07580 


Cxszx ITERATIVE SCHEME SXSBXSXSSXSSBBSSSSBBBSSBBBSSSSXBSBSSBSSXXSE 

BBBBBBSMS07590 


c 

SMS07600 


DO in ITsI ,MAXPAS 

SMS07610 


IAsI 

SMS07620 


SUMxO.O 

SMSO7630 


DO M2 J=3,NJM2 

SMS076* 4 


DO 42 I=3,NIM2 

SMS07650 


UIJxUIJO 

SMS07660 


C 

SMS07670 


CSSSB CHECK FOR OBSERVATION SSSSSSBSBSSSSBBSSSBSBESSBSBSSSSSBSSBBS 

bbbesbSMS07660 


C 

SMS07690 


IF (UO(I,J)) M3, 

SMS07700 


1(4 ALxO.O 

SMS07710 


GO TO 45 

SMS07720 


43 ALxAA 

SMS07730 


UIJsUIJ+AL 

SMS07740 


c 

SMS07750 


Css BB EQUATION FOR RESIDUAL SSSSSBSSBSSBSSBBSSSSSEBSBSSSBBSSSSSSBS 

bbbbbbSMS07760 


c 

SMS07770 


45 RESs-AL*UO(I. J)+UIJ*U(I I J)+UI1J1 I (U(X+1 ,J)+U(I~1 ,J)+U(I,J+1 )+U(I.JSMS07780 


1-1 ) )+ALF4B»(U(I-1 , J-1 )+U(I-1 , J+1 )*U(I+1 ,J-1 )+U(I*1 , J+1 ) ) 

SMSO7790 


RESxRES+ALF4«(U(I,J+2)+U(I,J-2)+U(I+2,J)+U(I-2,J)) 

SMS07800 


RLAXPsI ,/UIJ 

3MS07810 


c 

SMS07820 


Cssss CORRECT GUESS OF U sbsssbbssssssssssbsssbbbssbsssbsbssbbsbbbssbb*bSMS07830 


C 

SMS07840 


U(I,J)sU(I,J)-RLAXP»RES 

SMS07&50 


C 

SMS07860 


Cssss CHECK FOR APPROXIMATION SATISFIED AT ALL POINTS ssssssssssss 

SMS07870 


C 

SMSO78B0 


IF (DABS (RES) -ERR) 46,46,47 

SMS07890 


47 IAs2 

SMS07900 


46 SUMsSUM+RES»RES 

SMS07910 

• 

42 CONTINUE 

SMS07920 

* i 

STDsDSQRT(SUMZGRD) 

SMS07930 


GO TO (41,41), I A 

SMS07940 


41 CONTINUE 

SMS07950 


WRITE W,510)IT,STD 

SMS07960 


9 WRITE(7,51 1 ) IT 

SMS07970 


510 F0RMAT(1X,I3,E12.5) 

3MS07980 


511 FORMAT ( IX, 17HN0. OF ITERATIONS, 15) 

SMS07990 


500 FORMAT (8H WEIGHTS/5E15.2/28H STD DEVIATION OF RESIDUAL ) 

SMS08000 


RETURN 

SMS08010 


END 

SMS08020 
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Appendix B. Listing of SAS programs which carried out the estimation of 
yield and yield variances for corn and soybeans. The programs given 
are for the county level. Other levels were estimated in a similar 
manner. An example program for computation of covariances is also 
given. 
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* SAS PROGRAM FOR CORN YIELD ESTIMATION FOR EACH COUNTY IN IOWA. 

* BASED ON A USDA YIELD MODEL , WHICH USES LINEAR REGRESSION TECHNIQUES 

* AND METEOROLOGICAL PREDICTOR VARIABLES. 

* PROGRAM IS USED TO PREDICT 1978 YIELDS WITH A MODEL DEVELOPED USING 

* YIELD DATA FROM 1932 TO 1977, AND MET DATA FROM 1932 TO 1978. 

* WRITTEN BY CAROL JOBUSCH AT LARS, 1981. 

* \ 

DATA YLDMET; SET METCROP2. CTY ; 

DROP SACRES SPROD S YIELD; 

IF STRATUM « QQ 
CYLD * CYIELD; 

IF YEAR - 70 OR YEAR * 78 THEN CYLD ■ . ; 

TREND1-0; TREND2-0; 

IF YEAR > 40 THEN TRENDl -YEAR-40 ; 

IF YEAR > 60 THEN TRENDl-20; 

IF YEAR > 60 THEN TREND2 -YEAR-6 0 ; 

IF YEAR > 72 THEN TREND2-12; 

PCP TMP5-PCP5 *TMP5 ; 

PCP~TMP6-PCP6*TMP6 ; 

JUN_T_SQ=TMP6 *TMP6 ; 

JUL_P«PCP7 ; 

J UL_T_S Q* JUL_T_DT * JUL_T_DT ; 

LABEL-TRENDI-LINEAR TREND 1941-1960; 

LABEL TREND2 -LINEAR TRENC 1961-1972; 

LABEL PCP TMP5-MAY TEMP*PRECIP INTERACTION; 

LABEL PCP~TMP6 - JUNE TEMP*PRECIP INTERACTION; 

LABEL JUN~T SQ - JUNE TEMP DFN SQUARED; 

LABEL JUL_P = JULY PRECIPITATION DFN; 

LABEL JUL_T_DT - JULY TEMP DEPARTURE FROM TREND; 

LABEL JUL T_SQ = JULY TEMP DFT SQUARED; 

LABEL AUGJT_DT -AUGUST TEMP DEPARTURE FROM TREND; 

PROC REG DATA-YLDMET OUTSSCP-SSYX OUTEST-BDATA; 

TITLEl *******************************************************♦****. 
TITLE2 ******* IOWA CORN MODEL - COUNTY QQ ********************** 
TITLE3 ************************************************************* 
TITLE4 ; 

TITLES PREDICTION OF IOWA CORN YIELDS BASED ON 1932-1977 (EXCEPT 1970) ; 
MODEL CYLD-TRENDl TREND2 PCP_TMP5 PCP_TMP6 JUN T SQ JUL P JUL T_DT 
JUL_T_SQ AUG_T_DT/ P ; 

OUTPUT OUT-YLDMET PREDICTED-CPYIELD RESIDUAL-CRESID; 

PROC PRINT DATA-YLDMET; VAR YEAR CYIELD CPYIELD CRESID; 

PROC PLOT DATA-YLDMET; 

TITLES "'ACTUAL (*) VS PREDICTED (P) CORN YIELDS (1932-78)'; 

PLOT CYIELD* YEAR-'*' CPYIELD* YEAR- 'P' / OVERLAY ; 

* 

* PREPARE TO CALCULATE THE VARIANCE OF THE PREDICTED YIELD FOR 1978 

* ? 

DATA X78; SET YLDMET ; 

KEEP TRENDl TREND2 PCP_TMP5 PCP TMP6 JUN T SQ JUL P JUL T DT JUL T SQ 
AUG T_DT ; ~ ” 

IF YEAR-78; 
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******* 

****** 


DATA VSAVE ; SET YLDMET ; 

KEEP SYSTEM STRATUM CSIZE CYIELD CPYIELD CACHES CPROD; 

IP YEAR-78; 

DATA XPX; SET SSYX; 

IF _N_ -2 THEN DELETE; 

PROC MATRIX ; 

TITLE4 ******* 

TITLE5 ******* X"X MATRIX FOR THE 1978 ESTIMATE 

TITLE6 ******* *******. 

TITLE7 ************************************************************* 
FETCH XPX DATA- XPX ( KEEP- I NTERCEP TRENDl TREND2 PCPJTMP5 PCP TMP6 
JUN_T SQ JUL P JUL_T_DT JUL_T_SQ AUG_T_DT) COLNAME-XNAMES; 

FETCH X78“DATA»X78 (KEEP -TRENDl TREND2 PCPJTMP5 PCP_TMP6 JUN__T SQ 
JUL P JUL T DT JUL T SQ AUG T_DT) COLNAME-X7 8NAMES ; 

ONE-1 ; 

X78-ONE || X78; 

NAMEONE - "INTERCEP'* ; 

X78NAMES - NAMEONE | | X78NAMES ; 

PE-X78*INV(XPX) * (X78) “* ; 

FETCH SIGMA DATA-BDATA (KEEP*_SIGMA_) ; 

SIGMASQ * SIGMA*SIGMA ; 

VARCORN - SIGMASQ*PE ; 

IVARCORN ■ SIGMASQ* (1+PE) ; 

FETCH VSAVE 

DATA-VSAVE ( KEEP ‘SYSTEM STRATUM CSIZE CYIELD CPYIELD CACRES CPROD) 
COLNAME-VNAMES; 


VSAVE = VSAVE 
VSAVE - VSAVE 
NVARCORN 
VNAMES - 


VARCORN ; 
IVARCORN ; 

'VARCORN' 'IVARCORN' ; 
VNAMES | | NVARCORN ; 


PRINT XPX COLNAME-XNAMES ROWNAME-XNAMES ; 

PRINT X78 COLNAME-X7 8NAMES ; 

PRINT PE SIGMA ; 

PRINT VSAVE COLNAME-VNAMES; 

OUTPUT VSAVE OUT-SASOUT. CTYQQ COLNAME-VNAMES ; 

* 

* SAVE RESIDUALS FOR LATER CALCULATION OF THE COVARIANCE OF THE 

* PREDICTED YIELD FOR EACH STRATIFICATION SYSTEM. 

h • 

DATA CRESIDX. CTYQQ ; SET YLDMET; 

KEEP SYSTEM STRATUM YEAR CRESID CPYIELD ; 
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* SAS PROGRAM FOR SOYBEAN YIELD ESTIMATION FOR EACH COUNTY IN IOWA. 

* BASED ON A USDA YIELD MODEL, WHICH USES LINEAR REGRESSION TECHNIQUES 

* AND METEOROLOGICAL PREDICTOR VARIABLES. 

* PROGRAM IS USED TO PREDICT 1978 YIELDS WITH A MODEL DEVELOPED USING 

* YIELD DATA FROM 1932 TO 1977, AND MET DATA FROM 1931 TO 1978. 

* WRITTEN BY CAROL JOBUSCH AT LARS, 1981. 

* 

DATA YLDMET; SET METCROP2. CTY ; 

DROP CACRES CPROD CYIELD; 

IF STRATUM - QQ ; 

SYLD « S YIELD; 

IF YEAR - 78 THEN SYLD - . ; 

TREND* YEAR-31 ; 

IF YEAR > 74 THEN TREND*43; 

PCP_TMP 5 = PCP 5 * TMP 5 ; 

A UG _P_SQ* p CP8* p C p 8 ; 

LABEL TREND=LINEAR TREND 1932-1974; 

LABEL CUM_PCP=CUMULATIVE PRECIP OCT-APR DFN; 

LABEL PCP_TMP5=MAY TEMP*PRECIP INTERACTION; 

LABEL TMP 6 = JUNE TEMPERATURE DFN; 

LABEL PCP7 =JULY PRECIPITATION DFN; 

LABEL JUL_T_DT = JULY TEMP DEPARTURE FROM TREND; 

LABEL PCP 8 = AUGUST PRECIPITATION DFN; 

LABEL AUG__P__SQ * AUGUST PRECIPITATION DFN SQUARED; 

LABEL AUG_TJDT = AUGUST TEMP DEPARTURE FROM TREND; 

PROC REG DATA* YLDMET OUTSSCP=SSYX OUTEST=BDATA ; 

TITLE1 ************************************************************; 
TITLE2 ******* IOWA SOYBEAN MODEL - COUNTY QQ ******************; 
TITLE3 ************************************************************* 
TITLE4 ; 

TITLES BOOTSTRAP TEST FOR THE YEAR 1978; 

MODEL SYLD=TREND CUM_PCP PCP_TMP5 TMP 6 PCP7 JUL_T_DT PCP 8 AUG_P_SQ 
AUG_T_ DT ; 

OUTPUT OUT® YLDMET PREDICTED=SPYIELD RESIDUAL*SRESID; 

PROC PRINT DATA® YLDMET; VAR YEAR S YIELD SPYIELD ; 

TITLE5 PREDICTION OF IOWA SOYBEAN YIELDS BASED ON YEARS 1932-1977; 

PROC PLOT DATA® YLDMET; 

TITLE5 ''ACTUAL (*) VS PREDICTED (P) SOYBEAN YIELDS (1932-78)'; 

PLOT S YIELD* YEAR®'*' SPYIELD* YEAR* ' P ' / OVERLAY ; 

* 

* PREPARE TO CALCULATE VARIANCE OF THE PREDICTED YIELD FOR 1978. 

* 

DATA X78 ; SET YLDMET; 

KEEP TREND CUM_PCP PCP_TMP5 TMP 6 PCP7 JUL_T DT PCP8 AUG_P SQ AUG_T__DT; 
IF YEAR*78; 

DATA VSAVE; SET YLDMET; 

KEEP SYSTEM STRATUM S YIELD SPYIELD SACRES SPROD; 

IF YEAR=78 ; 

DATA XPX; SET SSYX; 

IF N_ *2 THEN DELETE; 
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PROC MATRIX ; 

TITLE4 ******* ******** 

TITLES ******* X""X MATRIX FOR THE 1978 ESTIMATE ******** 

TITLE6 ******* ******* ? 

TITLE7 ************************************************************* 

FETCH XPX DATA=XPX ( KEEP^INTERCEP TREND CUM_PCP PCP TMP5 TMP6 PCP7 # it 

JUL T DT PCP8 AUG P SQ AUG_T_DT) COLNAME=XNAMES ; 

FETCH x78”dATA«X78(KEEP» TREND CUM_PCP PCP_TMP5 TMP6 PCP7 JUL T DT 
PCP8 AUG P SQ AUG * DT) COLNAME«X78NAMES; 

ONE=l; 

X78=ONE || X78 ; 

NAMEONE * "INTERCEP" ; 

X78NAMES = NAMEONE | | X78NAMES ; 

PE»X78*INV(XPX) *(X78) 

FETCH SIGMA DATA=BDATA (KEEP=_SIGMA ) ; 

SIGMASQ = SIGMA*SIGMA ; 

VARSOY « PE*SIGMASQ ? 

IVARSOY * SIGMASQ* (1+PE) ? 

FETCH VSAVE 

DATA=VSAVE (KEEP=SYSTEM STRATUM SYIELD SPYIELD SACRES SPROD) 
COLNAME=VNAMES ; 

VSAVE = VSAVE VARSOY ; 

VSAVE - VSAVE IVARSOY ; 

NVARSOY * '‘VARSOY"' "IVARSOY"' ; 

VNAMES = VNAMES | | NVARSOY ; | 

PRINT XPX COLNAME=XNAMES ROWNAME=XNAMES f 
PRINT X78 COLNAME=X78NAMES; 

PRINT PE SIGMA } 

PRINT VSAVE COLNAME=VNAMES ; I 

OUTPUT VSAVE OUT=SOYOUT. ZZZ COLNAME=VNAMES ? I 

* 

* SAVE RESIDUALS FOR LATER CALCULATION OF THE COVARIANCE OF THE } 

* PREDICTED YIELD FOR EACH STRATIFICATION SYSTEM. 

* i i 

DATA SRESIDX. CTYQQ ; SET YLDMET; | 

KEEP SYSTEM STRATUM YEAR SRESID SPYIELD ; I 
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* SAS PROGRAM TO CALCULATE THE COVARIANCE OF THE YIELD ESTIMATION 

* FOR A GIVEN STRATIFICATION SYSTEM. 

* THE ACRES OF CORN AND SOYBEANS FOR EACH STRATUM ARE USED AS WEIGHTS. 

* WRITTEN BY CAROL JOBUSCH AT LARS, AUGUST 1981. 

* ? 

DATA WTEMPO ? SET METCROP2 . XXX ? 

IF YEAR « 78 ? 

PROC SUMMARY DATA-WTEMPO ? 

CLASS STRATUM? VAR CACRES SACRES; 

OUTPUT OUT-WTEMP1 SUM-CACRES SACRES? 

DATA WTEMP2? SET WTEMP1? RETAIN CTOT STOT? 

IF _TYPE -0 THEN DO? 

CTOT - CACRES? STOT - SACRES? DELETE? END? 

CWT - CACRES/CTOT? SWT * SACRES/STOT? SYSTEM * ZZ ? 

KEEP SYSTEM STRATUM CWT SWT? 

PROC SORT DATA-RESID. XXX OUT-TEMP? BY YEAR? 

DATA TEMP2 ? 

KEEP SYSTEM YEAR CRESIDl-CRESIDQQ 

SRESIDl-SRESIDQQ ? 

ARRAY CRESIDS (STRATUM) CRESIDl-CRESIDQQ ? 

ARRAY SRESIDS (STRATUM) SRESIDl-SRESIDQQ ? 

DO OVER CRESIDS? 

SET TEMP? BY YEAR? 

CRESIDS = CRESID ? 

SRESIDS * SRESID ? 

IF LAST. YEAR THEN RETURN? END? 

TITLE DATA SET TEMP? 

PROC CORR NOCORR COV OUT-CYTEMP (TYPE-COV) DATA* TEMP 2 ? 

VAR CRESIDl-CRESIDQQ ? 

TITLE DATA SET CVCORNY. XXX ? 

DATA CYTEMP2 (TYPE=COV) ? SET CYTEMP? 

IP TYPE = "COV' ? 

PROC TRANSPOSE DATA-WTEMP2 OUT-CORN WT PREFIX-CWT; VAR CWT? 

PROC MATRIX ? 

FETCH COVM DATA-CYTEMP2 (KEEP-CRESIDl-CRESIDQQ ) COLNAME-CNAMES ? 
FETCH CWT DATA-CORNWT (KEEP-CWTl-CWTQQ ) ? 

CWTD - DIAG(CWT) ? 

COVM = CWTD * COVM * CWTD? 

OUTPUT COVM OUT=CYTEMP3 COLNAME-CNAMES? 

DATA CYTEMP4? SET CYTEMP 3 ? 

SYSTEM * ZZ ? 

DATA CVCORNY. XXX ? SET CYTEMP 4 ? BY SYSTEM? 

KEEP SYSTEM COVCY? 

ARRAY CY(I) CRESIDl-CRESIDQQ ? 

I - 1? 

DO WHILE (I LT _N_) ? 

COVCY + CY? 

I + 1? 

END? 

IF LAST. SYSTEM THEN OUTPUT? 

PROC PRINT? 
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PROC CORR NOCORR COV OUT-SYTEMP (TYPE-COV) DATA=TEMP2 ; 

VAR SRESID1-SRESIDQQ ; 

TITLE DATA SET CVSOYY. XXX ; 

DATA SYTEMP2; SET SYTEMPj 
IF TYPE « "COV'* ; 

PROC TRANSPOSE DATA=WTEMP2 OUT=SOYWT PREFIX=SWT; VAR SWT; 

PROC MATRIX ; 

FETCH COVM DATA=SYTEMP2 (KEEP=SRESIDl-SRESIDQQ ) COLNAME*SNAMES ; 
FETCH SWT DATA=SOYWT ( KEEP=SWTl-SWTQQ ); 

SWTD = DIAG(SWT); 

COVM - SWTD * COVM * SWTD; 

OUTPUT COVM OUT=SYTEMP3 COLNAME=SNAM;3S ; 

DATA SYTEMP4; SET SYTEMP3; 

SYSTEM » ZZ ; 

DATA CVSOYY. XXX ; SET SYTEMP4 ; BY SYSTEM; 

KEEP SYSTEM COVSY; 

ARRAY SY(I) SRESIDl-SRESIDQQ ; 

1 * 1 ; 

DO WHILE (I LT _N_) ; 

COVSY + SY; 

I + 1; 

END; 

IF LAST. SYSTEM THEN OUTPUT; 

PROC PRINT; 
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Appendix C. FORTRAN programs used for estimation of the area variances. 
Both the pixel size (msefs3) and the field size (msefs) estimation 
programs are presented. 
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mse£s3 fortran 


IMPLICIT REAL* 8 (A-H,0~Z) 

REAL* 8 Y(100) ,YP(100) ,YPC(100) ,YW(100) , 

* E(100) ,YGPC(100) ,PE(100) , RE (100) ,Z(100) , 

* NEWPC,XXX(100) ,YYP(100) ,YAG(100) ,FLD(100) 

INTEGER* 4 JCNT (100) , STRATM , STRATO , COUNTY , SYSNUM , SNOLD 
REAL* 8 CTYNAM,CNAME (100) 

IEOF-Q 

1*0 

KKK*G 

READ (1,100, END=9 9 ) SYSTEM , S YSNUM , STRATM , COUNTY , CTYNAM , PLDACR , IGPC , 

* I PC , AGRI , CPACRE , CROP 
SYSOLD * SYSTEM 

SNOLD * SYSNUM 
STRATO = STRATM 
IP (IPC .NE. 0) GO TO 2 

1 READ (1,100, END“99) SYSTEM, SYSNUM, STRATM, COUNTY, CTYNAM, PLDACR, IGPC, 

* IPC, AGRI, CPACRE, CROP 

100 FORMAT (A8,2I3,I5,1X,A8,F4.0,2I3,2F8.0,2X,A8) 

IF ( IPC.EQ. 0) GO TO 1 

IF( (SYSNUM. NE. SNOLD) .OR. (STRATM. NE. STRATO) ) GO TO 3 !j 

2 1*1+1 l! 

CNAME ( I ) * CTYNAM n 

JCNT (I) “COUNTY f| 

Y ( I) =CPACRE d 

YAG ( I ) *CP ACRE /AGRI jj 

YPC (I) =IPC 

YGPC (I) =IGPC ,! 

FLD(I) = FLDACR II 

GO TO 1 If 

99 IEOF=l !j 

3 CONTINUE I 

KKK=KKK+1 If 

IF (I.,EQ. 1) GO TO 7 if 

SUM=0 . 0D00 If 

FLDSUM=0 . 0 f 

SUMPC=0 . 0D00 | 

SUMGPC=0 . 0D00 if 

DO 4 J=1,I if 

IF(YPC(J) ( ^?.n,0) SUM=SUM+Y(J) «! 

SUMPC=SUMFw+YPC ( J) 

SUMGPC*SUMGPC+YGPC ( J) I 

FLDSUM*FLDSUM+FLD ( J) 

4 CONTINUE 

IF (SUM. EQ. 0.0) GO TO 7 
P=SUM/ (SUMPC*25426 . 56D00) 

DO 55 J=1,I 1 

IF(YPC(J) .GT.0.0) YP(J)=Y(J)/(YPC(J) *25426. 56D00) i 

YW( J) =YPC ( J) /SUMPC }■ 

E ( J) * (P-YP ( J) ) **2 | 

55 CONTINUE i 
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FLDSIZ » FLDSUM/I 

A « P * (1.0-P) * 4.0 / 9.0 

NN»0 

DO 54 J*1,I 

IF(YPC(J) .EQ.0.0) GO TO 54 

NN-NN+1 

XXX (NN) *YPC ( J) 

YYP (NN) “E ( J) 

54 CONTINUE 

WRITE (2 ,225) CROP,SYSOLD,STRATO 
225 FORMAT ( '1 CROP ■ ' ,A8, 2X, 'SYSTEM »',A8,2X, ‘'STRATUM «',I3) 

DO 227 J«1,NN 

227 XXX (J)*XXX(J)* (25426. 0/FLDSIZ) 

CALL FITB (NN,XXX, YYP,A,B) 

A*A/( (FLDSIZ* (22932.0/25426.0) ) **B) 

SF»22932.0D00**B 
VAR ■ A*SF 

WRITE (2, 200) I,P,A,B,SF, VAR 

200 FORMAT (//' NUMBER OF COUNTIES IN STRATUM - ',13/' P «',F16.9 

* //' A ■' ,F16.9 /' B -',F16.9 

* /' SF »',F15.9 /' VAR »',F14.9) 

DO 11 J-1,I 

IF ( YPC ( J) . GT .0.0) PE(J) "A* (YPC(J) *22932. 0D00) **B 
IF(YPC(J) .EQ.0.0) PE(J) *0. 0D00 
IF (YPC ( J) .GT.0.0) RE(J)-(PE(J)-E(J) ) 

IF(YPC(J) .EQ.0.0) RE(J) «0 . 0D00 
11 CONTINUE 
M«0 

X-O.ODOO 
SX-0.0D00 
DO 13 J«1,I 

IF(YPC(J) .EQ.0.0) GO TO 13 

M«M+1 

X*X+RE(J) 

SX»SX+RE(J) **2 

13 CONTINUE 

IF (M.GT. 1) GO TO 14 
WRITE (2, 500) 

500 FORMAT (//, IX, 'DEGREES OF FREEDOM”0 ' , //) 

GO TO 7 

14 XM-M 
SD«SX-X**2/XM 
SD«SD/(XM-l.OD00) 

SD*DSQRT(SD) 

BAR-X/XM 

DO 17 J»1,I 

IF (YPC ( J) .GT.0.0) Z ( J) ■ (RE ( J) -BAR) /SD 
IF(YPC(J) .EQ.0.0) Z(J)«0.0D00 
17 CONTINUE 
WRITE (2, 400) 

400 FORMAT (//, IX, 'CNTY GPC PC ACRES PI WI', 

* ' (PI-P) **2 PROJECTED ERROR ', 

* ' Z VALUE CNTY') 


DO 6 J-1,I 

WRITE (2, 300) JCNT( J) ,YGPC(J) ,YPC(J) ,Y(J) ,YP(J) ,YW(J) ,E(J) , 

* PE(J) ,RE(J) , 1 ( 3 ) ,CNAME(J) 

300 FQRMAT(1X,13,2(3X,F4.0) ,3X,F8.0,2(3X,F7.4) ,2(3X»R10.7) , 

,1 (3X,F10.7) , (3X,F10.4) ,6X,A8) 

6 CONTINUE 

DO 21 J-1,I 
PC- YPC ( J) 

IF(YAG(J) .GE.1.0) GO TO 180 

IF ( (DABS(Z(J)) .LE.3.0) .OR* (PC.LE.0.0) .OR. 

* ( YPC ( J) . GT. 0 . 5*XGPC ( J) ) ) GO TO 21 
YPC ( J) *0.0 

WRITE(2,700) JCNT(J),Z(J) 

700 FORMAT ( // , IX, ''COUNTY'* , 15, 5X, "'REJECTED'* , 5X, 'Z«' ,F10. 3) 

GO TO 21 
180 CONTINUE 
YPC (J) “0.0 

WRITE (2, 780) JCNT(J) , YAG(J) 

780 FORMAT (// , IX , "COUNTY'' , 15 , 5X , ''CROP TO AG RATIO-' ,F10 . 3) 

21 CONTINUE 
NEWPC-0.0 
DO 22 J»1,I 

22 NEWPC»NEWTC+YPC(J) 

III-NEWPC 

JJJ-SUMPC 

WRITE (3, 1000) CROP , S YSOLD , STRATO , KKK , SUMGPC , SUMPC , 

* A,3,SF,P,VAR, I 

1000 FORMAT ( 2 A8 , 2 1 3 , 2F‘7 . 0 , 5F8 . 4 , 1 4 ) 

IF( (I.TX.LT. JJJ) .AND. (KKK.LE.3) ) GO TO 3 
WRITE (4,1000) CROP , S YSOLD, STRATO, KKK , SUMGPC, SUMPC 
6,A,B,SF,P, VAR, I 

7 CONTINUE 
SYSOLD-SYSTEM 
SNOLD ■ SYSNUM 
STRATO-STRATM 
1=0 

KKK-0 

IF ( IEOF.EQ. 0) GO TO 2 

STOP 

END 


SUBROUTINE PITB (N, X, Y, A, B) 

IMPLICIT REAL *8 (A-H,0-Z) 

REAL*8 X(100) ,Y(100) ,MAX 
EPS-0. 000001D00 
1-0 
K-10 
XK-K 
Bl— 0.90 
B2— 0.10 
WRITE ( 2 , 200) 

1 CONTINUE 

DELTA- (B2-B1)/XK 
MAX-P (X,Y,A,B1,N) 

DO 9 J-l , K 

B«B1+J*DELTA 

FB-F(X,Y,A,B,N) 

IF (MAX. LT. FB) GO' TO 9 

MAX-PB 

BB*B 

9 CONTINUE 
I-I+l 

WRITE ( 2 , 10 0 ) I , BB , MAX , DELTA 

IF(I.GT.20) GO TO 99 

I P ( DELTA . LT . EPS ) GO TO 99 

Bl-BB-DELTA 

B2-BB+DELTA 

GO TO 1 

200 FORMAT(//,T6,'K',T23,'B(K)',T41,'F(B()<))'rT64, 'DELTA') 
100 FORMAT (IX, 15,6 (1X,F20. 15) ) 

99 CONTINUE 
B-BB 
FB-MAX 
RETURN 
END 

REAL FUNCTION P*8 (X,Y,A,B,N) 

IMPLICIT REAL*8(A-H,0-Z) 

REAL* 8 X(100) ,Y(100) ,A,B,S,Y1,XB 

S-0.0D00 

DO 1 J-l ,N 

XB-1.0 

Yl-Y (J) 

IF( (X(J) .EQ.0.0) .OR. (B.EQ.0.0) ) GO TO 10 
XB-X ( J) **B 
10 S-S+ (Yl-A*XB) **2 
1 CONTINUE 
P-S/N 
RETURN 
END 


msefs fortran 


IMPLICIT REAL*8 (A-H,0-Z) 

GENERIC 

REAL*8 Y (100) ,YP(100) , YPC(IOO) , YW(100) , 

* E(100) , YGPC (100) ,PE(100) ,RE(100) ,Z(100) , 

* NEWPC,XXX(100) , YYP(IOO) ,YAG(100) ,FLD(100) 

INTEGER* 4 JCNT(IOO) , STRATM , STRATO , COUNTY , S YSNUM , SNOLD 
REAL* 8 CTYNAM,CNAME(100) 

IEOP*0 

1*0 

KKK*0 

READ (1,100, END-99) SYSTEM, S YSNUM, STRATM, COUNTY, CTYNAM, PLDACR, IGPC 

* IPC,AGRI,CPACRE,CROP 
SYSOLD * SYSTEM 

SNOLD « SYSNUM 
STRATO - STRATM 
IF (IPC .NE. 0) GO TO 2 

1 RE AD ( 1 , 1 0 0 , END- 9 9 ) S YSTEM , S YSNUM , STRATM , COUNT Y , CT YN AM , FLDACR , I GPC 

* IPC, AGRI,CPACRE, CROP 

100 FORMAT (A8,2I3,I5,lX,A8,F4.0,2I3,2F8.0,2X,A8) 

IF ( IPC.EQ. 0) GO TO 1 

IF ( (SYSNUM. NE. SNOLD) .OR. (STRATM. NE. STRATO)) GO TO 3 

2 I—I+l 

* CNAME(I) - CTYNAM 
JCNT ( I ) -COUNTY 
Y ( I ) -CPACRE 
YAG( I) -CP ACRE/AGRI 
YPC (I) -IPC 
YGPC (I) -IGPC 
FLD(I) - FLDACR 
GO TO 1 
99 IEOF-1 

3 CONTINUE 
KKK-KKK+1 

IF(I.EQ.l) GO TO 7 
SUM-0. 0D00 
FLDSUM-0 . 0 
SUMPC-0 . 0D00 
SUMGPC-0 . 0D00 
DO 4 J-1,I 

IF(YPC(J) .GT.0.0) SUM-SUM+Y ( J) 

SUMPC-SUMPC+YPC ( J) 

SUMGPC-SUMGPC+YGPC ( J) 

FLDSUM-FLDSUM+FLD ( J) 

4 CONTINUE 

IF (SUM. EQ. 0.0) GO TO 7 
P-SUM/ (SUMPC*25426 . 56D00) 

DO 55 J=1 , I 

IF (YPC ( J) . GT .0.0) YP ( J) =Y ( J) / ( YPC ( J) *25426 . 56D00) 

YW ( J) -YPC ( J) /SUMPC 
E ( J) * (P-YP ( J) ) **2 
55 CONTINUE 
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FLDSIZ - FLDSUK/I 
AA - P * (1.0-P) 

XP « FLDS 1 7, * 22932. / 25426,56 
NN*0 

DO 54 3*1,1 

IF(YPC(J) .EQ.0.0) GO TO 54 
NN-NN+1 

XXX (NN) *YPC ( J) * 22932. 

YYP (NN) *E ( J) 

54 CONTINUE 

WRITD(2* 225) CROP , S YSOLD , STRATO 
225 FORMAT ( ''I CROP * ' , A8 , 2X, "SYSTEM *',A8,2X, 'STRATUM -',13) 

A1 - P*(SQRT(XP)-1.2732)**2/XP 
A3 » P* (SQRT(XP) +1. 2732) **2/XP - A1 
A2 - 1. - A1 - A3 

A - Al* (l.-P) **2 + A2*P**2 -f A3*(0.3682-P+P**2) 

CALL FITB(NN,XXX,YYP,A,B) 

SF*22932. 0D00**B 
VAR ■ A*SF 

WRITE ( 2 , 200) I,P,A,B,SF,VAR 

200 FORMAT (//'* NUMBER OF COUNTIES IN STRATUM - ',13/' P »',F16.9 

* //' A «%F16.9 f* B ,F16, 9 

* /- SF # F15.9 /' VAR -',F14.9) 

DO 11 J«1,I 

IF ( YPC ( J) .GT.0.0) PE(J)«AMYPC(J) *22932 . 0D00) **B 
IF(YPC(J) .EQ.0.0) PE(J) “0.0D00 
IF (YPC ( J) .GT. 0.0) RE(J)-(PE(J)-E(J) ) 

IF(YPC(J) .EQ.0.0) RE ( J) «0 . 0D00 
11. CONTINUE 
M«0 

X-O.ODOO 
SX*0 . 0D00 
DO 13 J«1,I 

IF(YPC(J) .EQ.0.0) GO TO 13 

M*M+1 

X-X+RE(J) 

SX-SX+RE(J) **2 

13 CONTINUE 

IF (M.GT. 1) GO TO 14 
WRITE (2, 500) 

500 FORMAT (//, IX , 'DEGREES OF FREEDOM-O ',//) 

GO TO 7 

14 XM-M 
SD«SX-X**2/XM 
SD«SD/(XM-1.0D00) 

SD«DSQRT(SD) 

BAR-X/XM 

DO 17 J«1,I 

IF (YPC ( J) . GT .0.0) Z(J)»(RE(J)-BAR)/SD 
IF(YPC(J) .EQ.0.0) Z ( J) “0 . 0D00 
17 CONTINUE 


/ 
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WRITE (2 , 400) 

400 FORMAT (//, IX , 'CNTY GFC PC ACRES PI Wl', 

* ' (PI -P) **2 PROJECTED ERROR ', f 

* " Z VALUE CNTY') 

DO 6 J=1,I 

WRITE (2,300) JCNT ( J) ,YGPC(J) ,YPC(J) ,Y(J) ,YP(J) ,YW(J) ,E(J) , 

* PE(J) ,RE(J) ,Z(J) ,CNAME(J) 

300 FORMAT ( IX f 13 , 2 (3X,F4 . 0) , 3X,F8 . 0 , 2 (3X,F7 . 4) ,2 (3X,F10 . 7) , 

1(3X,F10.7) , (3X,F10.4) ,6X,A8) 

6 CONTINUE 

DO 21 J=1,I 
PC=YPC ( J) 

IF(YAG(J) .GE.1.0) GO TO 180 

IF ( (DABS ( Z ( J) ) .LE.3.0) .OR. (PC.LE.0.0) .OR. 

* (YPC (J) .GT. 0. 5*YGPC (J) ) ) GO TO 21 
YPC (J) =0 . 0 

WRITE (2, 700) JCNT(J),Z(J) I 

700 FORMAT (//, IX , ’'COUNTY'* , 15 , 5X, "REJECTED'* ,5X,''Z= B '',F10,3) 

GO TO 21 j 

180 CONTINUE J 

YPC (J) =0.0 

WRITE (2 ,780) JCNT(J) ,YAG(J) 

780 FORMAT (//, IX , '‘COUNTY’* , 15 , 5X, "CROP TO AG RATIO*' ,F10 . 3) | 

21 CONTINUE | 

NEWPC=0 . 0 

DO 22 0=1,1 !| 

22 NEWPC=NEWPC+YPC(J) J 

I1I=NEWPC 1 

JJJ=SUMPC 

WRITE(3,1000) CROP , S YSOLD , STRATO , SNOLD , SUMGPC , SUMPC , f 

* A,B,SF,P,VAR,I i 

1000 FORMAT (2A8, 213, 2F7.0,5F8. 4, 14) f 

IF( (III.LT.JJJ) .AND. (KKK.LE.3) ) GO TO 3 
WRITE (4, 1000) CROP, SYSOLD, STRATO, SNOLD, SUMGPC, SUMPC 
6 , A,B,SF, P, VAR, I ! 

7 CONTINUE 

SYSOLD=SYSTEM j 

SNOLD = SYSNUM \ 

STRATO=STRATM i 

1=0 \ 

KKK— 0 !i 

IF ( IEOF.EQ. 0) GO TO 2 | 

STOP I 

END I 


SUBROUTINE FITB (N, X, Y , A, B) 

IMPLICIT REAL* 8 (A-H,0-Z) 

REAL*8 X(100) ,Y(100) ,MAX 
EPS=0 .000001D00 
1 = 0 
K*10 
XK=K 
Bl=-0 .90 
B2=-0 . 10 
WRITE(2,200) 

1 CONTINUE 

DELTA® (B2-B1)/XK 
MAX=F(X,Y,A,B1,N) 

DO 9 J=1 , K 

B=B1+J*DELTA 

FB=F(X,Y,A,B,N) 

IF (MAX. LT. FB) GO TO 9 

MAX=FB 

BB=B 

9 CONTINUE 
1 = 1+1 

WRITE(2, 100) I,BB, MAX, DELTA 

IF ( I .GT. 20) GO TO 99 

IF (DELTA.LT. EPS) GO TO 99 

B1=BB-DELTA 

B2=BB+DELTA 

GO TO 1 

200 FORMAT(//,T6,'K',T23,'B(K) ',T41,'F(B (K))',T64, 'DELTA') 
100 FORMAT (lX f I5 r 6(lX r F20.15) ) 

99 CONTINUE 
B=BB 
FB=MAX 
RETURN 
END 

REAL FUNCTION F*8 (X,Y,A,B,N) 

IMPLICIT REAL* 8 (A-H r O-Z) 

REAL* 8 X(100) ,Y(100) ,A,B,S,Yl,XB 

S=0 . 0D00 

DO 1 J=1,N 

XB=1 . 0 

Yl=Y ( J) 

IF( (X(J) . EQ.0.0) .OR. (B.EQ.0.0) ) GO TO 10 
XB=X(J) **B 
10 S=S+ ( Y1-A*XB) **2 
1 CONTINUE 
F=S/N 
RETURN 
END 


